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Résumé : 



Ce rapport de recherche résume ravancement de travaux menés conjointement par l'IRCCyN et 
l'Ecole Polytechnique de Montréal concernant la résolution du problème inverse pour l'imagerie 
sismique des fondations de pylônes électriques. Nous abordons plusieurs méthodes de type « carto- 
graphie ». Nous nous intéressons plus particulièrement à des méthodes basées sur une formulation 
bilinéaire du problème direct d'une part (CSI, gradient modifié, etc.) et à des méthodes basées sur 
une formulation dite « primale » d'autre part. Les performances de ces méthodes sont évaluées 
sur des données synthétiques. 

Ces travaux ont été partiellement financés par RTE - CNER, qui est à l'initiative du projet, et 
ont été effectués avec la collaboration de EDF R&D. 
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Introduction 



Le problème d'imagerie des fondations de pylônes consiste à déterminer la forme d'objets enfouis dans le sol. Pour 
cela, on génère des ondes sismiques à l'aide d'une source placée en surface et on mesure une partie du champ de vitesse 
résultant à l'aide d'une série de capteurs également placés en surface. Plusieurs tirs sont effectués en modifiant la position 
de la source d'un tir à l'autre. 

Dans un premier temps, nos deux équipes (l'IRCCyN à Nantes et l'Ecole Polytechnique de Montréal) ont abordé ce 
problème en travaillant de façon commune sur une première méthode d'inversion (la méthode CSI). Ce travail a fait 
l'objet du précédent rapport d'avancement [ ]. Nous nous sommes ensuite partagé les tâches de manière à travailler sur 
deux familles de méthodes différentes. Le tableau ci-dessous détaille les points abordés par chacune des équipes. Afin 
de pouvoir comparer les performances des différentes méthodes abordées, nous les avons testées sur les mêmes jeux de 
données synthétiques. 



Jusqu 'en juillet 2009 


- Travail commun entre Nantes et Montréal - 

• Prise en main de l'algorithme de résolution du problème direct 

• Recherche bibliographique sur la méthode CSI 

• Adaptation de la méthode CSI à notre problème 

• Début de l'implémentation de la méthode 


Jusqu'en septembre 2009 


- Partage des tâches entre Nantes et Montréal - 

A l'Ecole Polytechnique de Montréal : 

• Poursuite du travail concernant la méthode CSI 

• Proposition de plusieurs variantes 
A l'IRCCyN : 

• Recherche bibliographique sur les méthodes adaptées à la formulation "primale" 


Jusqu'en décembre 2009 


A l'Ecole Polytechnique de Montréal : 

• Travail sur la formulation du problème direct 

• Implémentation des variantes de la méthode CSI et premiers résultats 

• Travail sur le problème de conditionnement 
A l'IRCCyN : 

• Travail sur l'optimisation de la formulation primale 

• Implémentation des méthodes de type gradient et premiers résultats 

• Utilisation du changement de variables et premiers résultats 



Nous avons conservé les hypothèses formulées dans le rapport d'avancement précédent [l]. On rappelle notamment les 
points suivants : 

- l'approche abordée est de type "cartographique". Nous cherchons donc à reconstruire des cartes du sous-sol sans 
utiliser d'information a priori sur la forme de l'objet recherché ; 

- nous considérons que la distribution de deux caractéristiques mécaniques (les vitesses des ondes de pression et 
cisaillement) est suffisante pour décrire un milieu ; 

- la méthode de résolution du problème direct, qui permet de construire un jeu de données synthétiques pour un milieu 
donné, est basée sur les différences finies. 

Dans la première partie de ce rapport, nous reviendrons sur la formulation du problème direct. Une étude approfondie 
nous a permis de mettre en évidence une structure particulière de la matrice dite "d'impédance" faisant intervenir des 
opérateurs de filtrage. 

La seconde partie concerne les différentes méthodes d'inversion. Nous reviendrons dans un premier temps sur la méthode 
CSI puis nous présenterons les autres méthodes d'inversion que nous avons abordées. Les premières sont analogues à la 
méthodes CSI et utilisent une formulation faisant intervenir des variables auxiliaires. Les autres se basent sur la formulation 
dite "primale", il s'agit de méthodes n'agissant que sur les variables d'intérêt (pas de recours à des variables auxiliaires). 
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Chapitre 1 

Le problème direct 



1.1 Les équations de propagation 



Le problème direct consiste à construire des données synthétiques en modélisant la propagation des ondes dans un 
milieu donné. On utilise pour cela les équations de propagation en milieu élastique. Ces équations sont présentées dans 
le domaine fréquentiel (on note u) la pulsation de la fonction excitatrice et Nf le nombre de fréquences considérées) et 
sont discrétisées à l'aide des différences finies, ce qui permet d'écrire le problème sous la forme d'un système d'équations 
linéaire (voir [2] et [j] pour plus de détails). Afin d'alléger les notations, nous considérerons dans un premier temps qu'il 
n'y a qu'une seule position de tir et que nous travaillons en monofréquentiel. 

Nous schématisons sur la Figure f .f le domaine d'étude. Le milieu de propagation D de dimensions finies est entouré 
d'une zone "PML", ce qui permet de simuler la propagation des ondes dans un milieu aux dimiension infinies (les PML 
atténuent les ondes de manière à éviter une réfiexion sur les bords). 




□ Milieu D 

■ Source 

■ Capteurs 

■ PML 

■ Bords 



Figure 1.1 - Description du modèle proposé 



A l'intérieur du milieu D et dans la zone PML (voir Figure 1.1), les équations de propagation dans le domaine fréquentiel 
s'écrivent : 
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Chapitre 1. Le problème direct 



-iujj.j,{r) = uj'^Vx{r) + arc{r,uj)dx{vp{rYa:c{r,uj)dxVj,{r)) 

+ a^{r,uj)d^{{vp{rY -2v,{rY)ay{r,uj)dyVy{r)) 
+ ay{r,uj)dy{vs{r)'^ax{r,uj)dxVy{r)) 
+ ay{r,u:)dy{vs{rfay{r,uj)dyV^{r)) 

-iujfy{r) = uj^Vy{r) + ax{r,uj)dx{vs{r)'^a^{r,uj)d^Vy{r) 
+ a^{r,u!)dxivs{r)'^ay{r,u!)dyVa;{r)) 
+ ay{r,uj)dy{{{vp{rf - 2v,{rf)a,{r,uj)d,V^{r)) 
+ ay{r,u})dy{{vp{rfay{r,uj)dyVy{r)) 

où fx{f) et fy{r) sont les composantes de la fonction source au point r (rapport de la force sur la masse volumique), Vx{r) 
et Vy{r) sont les composantes de la vitesse au point r, Vp{r) et Vs{r) sont les caractéristiques recherchées (respectivement 
les vitesses des ondes P et S) et ax{r,uj) et ay{r,uj) sont des coefficients introduits par les PML qui sont égaux à 1 à 
l'intérieur du domaine. 

Sur les bords du domaine, on a les contraintes suivantes : 

Vx{r)^0 et Vy{r) = (1.2) 

1.2 Discrétisation des équations de propagation 

1.2.1 Le stencil de Saenger 

Ces équations de propagation sont discrétisées en utilisant le stencil proposé par Saenger [■'!] et en introduisant des 
opérateurs de différences finies adaptés. On utilise deux grilles de points de résolution Ax suivant l'axe x et de Ay suivant 
l'axe y disposées en quinconce et définies sur l'ensemble du milieu. 




\ 

i=N,j=M 

Figure 1.2 - Grilles utilisées pour la discrétisation : la grille bleue est associée aux champs de forces {fx et fy) et de 
vitesses (T4 et Vy) et la rouge aux champs des caractéristiques {vp et Vs). Les coefficients associés aux PML {ax et ay) 
sont définis sur les deux grilles, ils sont différents de 1 uniquement dans la zone PML (surface grisée). 



La première grille (représentée par les points bleus sur la figure 1.2) est fiée aux champs de force fx et fy (source) et 
de vitesse T4 et Vy. Les points de cette grille, de coordonnées Xi = iAx et yj — jAy avec i et j entiers, sont désignés par 
les indices i et j. 

La seconde grille (représentée par les points rouges sur la figure 1.2) est liée aux caractéristiques du milieu Vp et Vg- Les 
points de cette grille, de coordonnées Xi — {i ± |)Ax et yj = (j ± |)Ay avec i et j entiers, sont désignés par les indices 
î ± i et j ± i . 
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Chapitre 1. Le problème direct 



Les coefficients ax{r,uj) et ay(r, w) sont quant à eux définis sur les deux grilles. 

Remarque : Si la première grille est constituée de M lignes et de N colonnes (si l'on reprend les notations de la 
Figure 1.1, on a A/ = Mt + 2Mpml + 2 et N = Nt + 2Npml + 2), la seconde grille est constituée de M — 1 lignes et de 
— 1 colonnes. 



1.2.2 Introduction des opérateurs de différences finies 

Nous allons reprendre les équations différentielles de propagation (Equations 1.1 et 1.2) et y introduire quatre opérateurs 
de différences finies adaptés au stencil de Saenger correspondant aux dérivées partielles spatiales dx et dy : 
- un opérateur Q'^ correspondant à la dérivée partielle selon x d'une grandeur associée à la première grille : 



1 



[(A)i+ij+i + (A)i+ij - (A)jj+i - iX)ij] (1.3) 



'^+2.3+2 2Ax 

un opérateur correspondant à la dérivée partielle selon y d'une grandeur associée à la première grille 



0iti;2j+i;2 
+ 



+1j 



1 



(A)i+ij + (x)ij+i - (a:)j,. 



(1.4) 



'^+2^+2 2Ay 

un opérateur "H^ correspondant à la dérivée partielle selon x d'une grandeur associée à la seconde grille 



^iti;2j+i;2 



|tti;2j-i;2 



1 



2Ax 



(X), 



{X\-l 



{X),_.,. (1.5) 



i-1/2,j+1/2 ^ ^i+1/2,j+1/2 



un opérateur correspondant à la dérivée partielle selon y d'une grandeur associée à la seconde grille 



(1.6) 



i-1/2J-1/2 ^ ^K1/2.j-1/2 
• l.] 

i-1/2j+1/2 ^ ^i+i;2jt1/2 



Remarque : L'utilisation de l'un de ces quatre opérateurs induit le passage d'une grille à l'autre. 
A l'intérieur du milieu D et dans la zone PML, on obtient alors les équations suivantes : 



(F,),,, = uo^{Vx\,j + (a.H"(îi>.a"(T4))).,j 

+ {ax'H-{{vl-2vl)aygy{Vy))),^, 
+ {ayHy{vlaxQ%Vy))),,, 

+ {ayHy{vlaygy{Vx)))^,, 
+ {axn^{v^saygy{Vx))) 



+ {ayny{{vl-2vl)axg''iVxm,, 



Sur les bords du domaine, on impose : 



+ {ayW{v'^ayQy{Vy))\^ 



et {Vy\^ = 



(1.7) 



(1.8) 



1.3 Ecriture sous la forme d'un système linéaire d'équations 

A partir de l'ensemble des équations discrétisées (Eq. 1.7 et 1.8), on forme un système linéaire d'équations reliant les 
forces aux vitesses que l'on écrit de la façon suivante : 



F = K .V 



(1.9) 



D. Vautrin, M. Voorons, J. Hier, Y. Goussard 
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ou, en distinguant les composantes horizontales et verticales des champs 







A 


A 




'Vx 






A 

^yx 


A 

^yy_ 




yy. 



(1.10) 

Si la première grille est constituée de M lignes et de N colonnes, les vecteurs F et V sont de longueur 2MN et la 
matrice A, appelée "matrice d'impédance", est une matrice bande de taille 2MN x 2MN. Ses coefficients sont déduits de 
l'ensemble des équations discrétisées ; ils dépendent donc des caractéristiques du milieu de propagation (champs Vp et v^) 
ainsi que des coefficients liés aux PML {ax{r,Lo) et ay{r,uj)). On donne ci-dessous les expressions des quatre sous- matrices 
qui la constituent : 



A,, = [A^,, + Diag{aï}H^Diag{î;|}Diag{a^}G'' + Diag{a?}HyDiag{î;2}Diag{a^}G''] 
A,y = [Diag{aî}H^Diag{î;2 _ 2î;2}Diag{a|}Gy + Diag{anH''Diag{î;2}Diag{a^}G^] 
Ay^ = [Diag{aîf}H^Diag{î;2}Diag{a|}Gy +Diag{af}HyDiag{î;^ - 2î;2}Diag{af }G^] 
Diag{aî}H''Diag{î;2}Diag{a^}G'' + Diag{af }HyDiag{î;^}Diag{a|}Gy] 



^yy ~ [-'^".a 



où Diaglw} {w étant un vecteur) désigne une matrice diagonale dont la diagonale est le vecteur w. 
Les matrices A^^^^ et A^^ ,^ sont identiques. Ce sont des matrices diagonales telles que : 



A^Aj + iM-l)t;j + iM-l)t) 



,(j + (A/-l)ï;j-H(M~l)î) 



1 si (i,j) appartient aux bords 
cj^ sinon 



(1.11) 
(1.12) 
(1.13) 
(1.14) 



(1.15) 



Les matrices G'', G^, H'^ et sont des matrices rectangulaires associées aux opérateurs de différences finies. 

- G^ et G^ sont de taille (M — 1){N — 1) x MN. Chaque ligne de ces deux matrices contient les coefficients de 
l'opérateur de différences finies associé (i^Sï -'-2Sï/)' 

- H'^ et H7 sont de taille MN x (M— 1)(A^— 1). Certaines de leurs lignes ne contiennent que des zéros (lignes associés 
aux points du bord du domaine) ; les autres contiennent chacune les coefficients de l'opérateur de différences finies 
associé (±2Sï et ±2^). 

Afin de prendre en compte les différentes pulsations qui constituent la fonction excitatrice et les différentes positions 
de tir de la source (Nk tirs), nous écrirons maintenant le système linéaire reliant le champ des forces au champ des vitesses 
de la façon suivante : 

(1.16) 

7^Nf) et k désigne à la position de la source 



^1, 



où uj correspond la pulsation de la fonction excitatrice considérée {ui 
ik=l,...,Nk). 

La résolution du problème direct consiste à résoudre le système linéaire donné par l'équation I.IG pour chaque fréquence 
et chaque position de la source et à retenir les composantes du champ de vitesse mesurées par les capteurs. 



1.4 Décomposition de la matrice d'impédance sous forme d'une somme 

La matrice d'impédance A^^^p^s du système linéaire 1.16 peut être décomposée en la somme de trois termes : 

A^,p,, = A^ + AP + A^ (1.17) 

où : 



- A^* contient les éléments de A^ p g qui dépendent des {vp)ij ; 

- A'' contient les éléments de A^^^p s qui dépendent des (ws)ij ; 

- Ac^ contient les éléments qui ne dépend ni des {vp)ij, ni des {vs)i,j (cette matrice ne dépend que de la fréquence w). 
En reprenant les expressions données dans la partie précédente, on obtient les expressions suivantes : 



AP = 



A" = 



Diag{af}H^' 
Diag{af}Hy 

DiagjaflHy 
Diag{aif}H^ 



-DiagjaflHy 



Diag{î;2} [Diag{af}G'' Diag{a|}Gy] 
Diag{î;2} [Diag{a^}Gy Diag{ai}G''] 
Bia.g{vl} [Diag{af}G- -Diag{a|}Gy] 



(1.18) 



(1.19) 
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-Diag{<}H''' 
-Diagjaj'IHy 



Diag{î;f} [Diag{af}G'' Diag{a^}Gy] 



La matrice A^^ est diagonale ; on a 



A, 







(1.20) 



Cette décomposition sera utilisée lors de la construction des équations utilisées pour résoudre le problème inverse. 

AP et A'' sont des matrices bandes. Si l'on construit les vecteurs F^^ ]^ et V^^k en alternant composantes horizontales 
et composantes verticales, la strucure des matrices A^ et A* est telle que représentée sur la Figure 1.3. On remarque que 
chaque ligne de ces matrices comporte ou 18 coefficients non nuls. En effet, à l'intérieur du domaine D et dans la zone 
PML, chaque composante du vecteur F^^k en un point du milieu s'écrit en fonction des deux composantes de la vitesse au 
même point et à ses huit voisins. Il s'agit en fait d'un voisinage du deuxième ordre puisque l'expression d'une composante 
de F^^)^ en un point passe par l'utilisation des opérateurs C, G^, et qui induisent le passage d'une grille à l'autre 
par une combinaison linéaire en quatre points voisins (cf équations 1.3 à 1.6). Sur les bords, la relation entre Fi^_k et V^i^k 
est donnée par la matrice A^^. 



2xM 



2xMxN 





2xM 
2 

2(M-2) 



2(M-2) 



2(M-2) 
2 

2xM 



2xMxN 

Figure 1.3 - Structure des matrices A'^ et A* si l'on construit F^ j^ et ^ en alternant composantes horizontales et 
composantes verticales. Chaque ligne comporte ou 18 coefficients non nuls. 



Remarque : Les coefficients des matrices A^* et A* associés à la zone PML dépendent des ax{(S) et ay{uS). Une partie 
de ces deux matrices est donc dépendante de la fréquence considérée uj. 
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Chapitre 2 



Introduction sur les méthodes de 
résolution du problème inverse 



Nous abordons maintenant le problème inverse. L'objectif est de retrouver la distribution des grandeurs caractéristiques 
d'un milieu en utilisant les données mesurées par des capteurs et la fonction source. Pour réduire la sous détermination 
du problème, nous utilisons les données issues de différentes positions de tir. Nous nous limiterons à la recherche de deux 
caractéristiques (vp et Vs) et nous considérerons que les données mesurées correspondent à la vitesse verticale Vy au niveau 
des capteurs situés en surface. 

Nous proposons plusieurs méthodes d'inversion qui ont toutes en commun le fait de minimiser un critère de façon 
itérative. On peut regrouper ces différentes méthodes en deux familles : 

- la première comprend les méthodes utilisant une formulation faisant intervenir des variables auxiliaires. L'avantage 
de ces formulations est qu'elles permettent de se ramener à la résolution de plusieurs sous-problèmes plus simples 
(minimisation d'un critère quadratique par exemple). Cependant, ce type de formulation nécessite l'optimisation 
d'un plus grand nombre de variables (variables d'intérêt et variables auxiliaires). 

- la seconde est associée aux méthodes n'utilisant pas de variable auxiliaire. Ces méthodes ont l'avantage de n'agir que 
sur les variables d'intérêt et donc de travailler dans un espace de dimension plus réduite. Cependant, la formulation 
utilisée est non linéaire ce qui rend le processus d'optimisation plus complexe. 

Pour la plupart des méthodes proposées, la recherche des caractéristiques physiques du milieu consistera à retrouver 
leurs variations par rapport à un milieu de référence. Le milieu de référence est choisi par l'utilisateur ; il s'agit, par exemple, 
d'un miUeu composé uniquement du matériau qui entoure l'objet difîractant (dans notre étude, les caractéristiques du 
milieu de référence sont celles de la terre). 

Les formulations établies ici découlent de l'équation matricielle associée au problème direct et de la propriété de 
décomposition de la matrice d'impédance sous forme de somme (voir Partie L4 page 12). 

2.1 Construction d'une première formulation bilinéaire 

2.1.1 Expression des vitesses en fonction des caractéristiques du milieu recherché 

Plaçons nous tout d'abord dans le milieu de référence qui est choisi par l'utilisateur. Pour désigner ses caractéristiques, 
nous utiliserons les notations Vp^ et Vs,o ■ Pour une fréquence et une position de la source données, le champ de vitesse 
qui se propage dans ce milieu de référence est dit "incident" ; on le note j,. On a la relation : 



où la matrice d'impédance {Ai^_p^s)o est liée aux caractéristiques du milieu de référence Vp^o et Vs,o (cette matrice est donc 
connue) ; elle peut être décomposée en la somme de trois termes comme expliqué dans la Partie 1.4. 

Considérons ensuite le milieu recherché, c'est-à-dire le milieu en présence de l'objet diffractant. Il est caractérisé par 
les champs Vp et v^. Dans ce milieu et pour une fréquence et une position de la source données, le champ de vitesse est 
dit "total" ; on le note Vuj.k- 



[A^+Ag + Ag].F^% 



(2.1) 



-A. 




(2.2) 
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Chapitre 2. Introduction sur les méthodes de résolution du problème inverse 



En soustrayant (2.1) à (2.2), on obtient : 

Vu,k = - {A^,p,s)ô\[A^ ~ Ag].K.,, + [A^ - Ag].K.,,) (2.3) 

Nous avons donc établi une relation entre le champ des vitesses Vuj^k et les caractéristiques du milieu recherché qui 
interviennent dans les matrices et A''. 



2.1.2 Introduction des contrastes et des sources de contraste 

Dans cette première formulation, nous nous intéressons aux variations des caractéristiques du milieu recherché par rap- 
port au milieu de référence. Nous chercherons plus particulièrement à retrouver la distribution spatiale de deux contrastes : 

- {Xp)i,j — (î'n — ^p,o)ij 6st le contraste en vitesse des ondes P ; 

- {Xs)i,j = {Vg ~ o)ij ^st le contraste en vitesse des ondes S. 

Ils caractérisent le milieu en présence de l'objet diffractant au même titre que les champs Vp et Vg puisque le milieu de 
référence est connu. Xp et Xs sont des vecteurs de longueur {M — 1){N — 1) (voir Partie 1.2.1). 

Les deux contrastes recherchés interviennent dans l'équation 2.3 via les matrices A^ — Aq et A* — Ag. Cela nous amène 
à définir deux matrices de contraste : 

- Xp = A^" — Aq est la matrice de contraste en vitesse des ondes P ; chacun de ses éléments est une combinaison 
linéaire des {xp)i,j ; 

- Xg = A" — Aq est la matrice de contraste en vitesse des ondes S ; chacun de ses éléments est une combinaison linéaire 
des {Xs)i,j- 

Xp et Xs sont des matrices de taille 2MN x 2MN. En reprenant les expressions données dans la partie 1.4, on obtient 
les expressions des matrices de contraste en fonction des contrastes : 

Diag{xp} [Diag{af IG'' Diag{af}Gy] (2.4) 



XP = 



Diag{aî}H^ 
Diag{af}Hy 



Diag{af}Hy 
Diag{aî}H'' 

Diag{af}H'' 
-DiaglafjHy 

-Diag{<}H^ 
-DiagjckflHy 



Diag{x.} [Diag{a|}Gy Diag{ai}G-] 
Diag{xJ [Diag{ai}G- -Diag{a|}Gy] 
Diag{xJ [Diag{ai}G'' Diagla^jG^] 



(2.5) 



On introduit finalement les sources de contraste : W^^^k — O^p + ^s) H s'agit des variables auxiliaires introduites 
pour construire une première formulation bilinéaire. Ce sont des vecteurs de taille 2MN. 
L'équation (2.3) peut alors s'écrire de la façon suivante : 



(2.6) 



2.1.3 Construction des équations de données 

Les équations de données correspondent à l'expression des données mesurées au niveau des capteurs en fonction 
des variables auxiliaires, pour chaque fréquence et chaque position de la source. Pour établir ces équations, on utilise 
l'expression de Vui^k précédente (Eq. 2.6) en ne retenant que les composantes mesurées par les capteurs : 

(Kj,fe)capt = (KÎj,fe)capt — Ei(A^^p^s)q FFi^^fc (2.7) 

OÙ : 

- El est une matrice d'échantillonnage de taille x 2MN. Elle correspond à l'ensemble des Ne lignes de la matrice 
identité associées aux composantes mesurées par les capteurs. 

~ ( Kj,fc)capt = El Vi^^k (composantes du champ de vitesse total mesurées par les capteurs) 

- ( fc)capt = El j, (composantes du champ de vitesse incident mesurées par les capteurs) 
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2.1.4 Construction des équations de couplage 



Les équations de couplage correspondent à l'expression des variables auxiliaires (les sources de contraste W^j^k) en 
fonction des variables d'intérêt (les contrastes Xp et Xs)- Pour établir ces équations, on reprend l'équation (2.6) et on 
multiplie les termes de droite et de gauche par la matrice Xj, + Xs . On détermine ainsi une équation de couplage pour 
chaque fréquence et chaque position de la source : 



(Xp + X,)(l/«,,-(A^,p,,)oiW^^,fc) 



(2.8) 



2.1.5 Réduction du problème à une zone d'étude 

Une réduction de la taille du problème est possible si l'on s'intéresse uniquement à un sous-domaine du milieu D 
(domaine situé à l'intérieur des PML, voir Figure 1.1 page 9). On définit une zone appelée "zone d'étude" à l'intérieur 
de laquelle on autorise les paramètres caractéristiques {vp et Vg) à prendre des valeurs différentes de celles du milieu de 
référence. Autrement dit, la zone d'étude est la région à l'intérieur de laquelle les contrastes Xp et Xs sont susceptibles de 
prendre des valeurs non nulles ; à l'extérieur de cette zone, les contrastes restent nuls. 

Il en résulte que plusieurs lignes et colonnes des matrices Xp et X^ sont nulles. On peut donc se ramener à deux 
matrices de taille réduite de taille iVo x A'o, où Nq désigne le nombre de lignes et de colonnes non nulles de Xp et Xj : 



(Xp)i.od — E2XpE2 (et donc Xp — E2(Xp)i.cdE2) 



(X.) 



rcd 



EaXsE* (et donc X^ = E* (Xs)rcdE2) 



(2.9) 
(2.10) 



oîi la matrice d'échantillonnage E2 est utilisée pour éliminer les lignes et les colonnes nulles des matrices de contraste. Elle 
est de taille A'o x 2MN et correspond à l'ensemble des iVo lignes de la matrice identité associées aux lignes non nulles des 
matrices Xp et X^. 

Étant donné que W^j^k — (Xp + X,,) Vi^^k, certains coefficients de Wu:,k sont nuls. On peut donc éliminer les coefficients 
nuls et se ramener à des vecteurs de source de contraste de taille réduite : 



(W^a;,fe)red = E2 W^^^fc (et douC W^^k = El{W^,k)rcd) 



(2.11) 



On peut alors réécrire les équations de données et de couplage en faisant intervenir les matrices de contraste et les 
vecteurs de source de contraste de taille réduite : 



(Kj,fc)capt — ( K!!,fe)capt — Ei(Atj^p^s)o ^E2( IVcj,fc)r 
( W^w,fc)rcd 



((Xp),.ed + (X,)„d) (E2 - E2(A^,p,,)o iE^( W^^,fe) 



rcd j 



(2.12) 
(2.13) 



La réduction de la taille du problème a également un impact sur l'expression des matrices de contraste en fonction des 
contrastes ((xp)zE et (xs)ze désignent les composantes des vecteurs de contraste appartenant à la zone d'étude) : 



(Xp) 



■p Jrcd 



Diag{(xp)zE} [G'' G^] 



(2.14) 



(X.) 



rcd 



Diag{(xs)zE} G''] 



Diag{(x.)zE} [G- -G^] 



-m 



Diag{(x.)zE} [G- G^] (2.15) 



011 les matrices G^, G'', H'' et H*' sont de taille réduite. 

Remarque : Les expressions ci-dessus ne font plus apparaître les coefficients af, a\, a| et a\ car la zone d'étude se 
situe à l'intérieur du milieu D (pas d'intersection avec la zone PML) ; ces coefficients sont donc égaux à 1. 



2.1.6 Simplification des écritures 

Afin de simplifier les écritures, et en tenant compte des remarques précédentes concernant la réduction des matrices 
et des vecteurs, nous écrivons désormais les équations de couplage et de données de la façon suivante : 



(T4;,/c)capt — ( Kj,fc)capt — Bf, FFtj^fc 

W^^k = (Xp + X,)(l/Ofe-B^IF^,,) 



(2.16) 
(2.17) 



où = E2(A^,p,,)o lE* et = Ei(A^,p,,)o ^E*. 
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Nous avons donc construit une formulation bilinéaire où les équations de données donnent l'expression des données 
synthétiques (( Vl;,fc)capt) eu fonction d'un jeu de variables auxiliaires (W^j^k) et où les équations de couplage font le lien 
entre les variables auxiliaires et les variables d'intérêt (xp et Xs)- 

Nous simplifions également l'expression des matrices de contraste en fonction des contrastes : 

Xp = HPDiag{xp}GP (2.18) 

3 

X, - ^H?Diag{x.}G? (2.19) 

2.2 Construction de la formulation primale 

La seconde famille de méthodes abordée correspond aux méthodes utilisant la formulation dite "primale", c'est-à-dire 
l'expression des données synthétiques en fonction des variables d'intérêt sans faire intervenir les variables auxiliaires. Cette 
relation peut se déduire facilement de l'expression bilinéaire proposée précédemment. 

Isolons tout d'abord le terme Wui^k dans l'équation de couplage (Eq. 2.17), pour une fréquence w et une position k de 
la source données : 

= [/+(Xp + X,)B^]-i(Xp + X,)F:;^, (2.20) 
Reportons cette expression dans l'équation de données (Eq. 2.16) : 

( K,,fe)capt = ( l^°,fc)capt - B^[/ + (Xp + X,)B^]-i(Xp + X,) V^^, (2.21) 

Nous obtenons ainsi l'expression des composantes de la vitesse mesurées par les capteurs ( Vtj,fc)capt en fonction des 
contrastes recherchés (xp et Xs)- 



2.3 Procédures de minimisât ion utilisées 

Pour les différentes méthodes d'inversion envisagées, l'objectif est de minimiser un critère C que l'on exprime en fonction 
des contrastes Xp et Xs ainsi que plusieurs variables auxiliaires pour certaines des méthodes abordées. Nous utilisons pour 
cela une méthode locale de type gradient qui consiste en une succession de minimisations en une dimension : à chaque 
itération, on définit une direction de recherche dans l'espace de représentation puis on détermine un pas de progression 
efficace le long de cette direction de recherche. 

2.3.1 La définition de la direction de recherche 

Dans le cas où le critère à minimiser est quadratique, (cas rencontré pour certaines méthodes utilisant des va- 
riables auxiliaires) nous utilisons le gradient conjugué linéaire. Pour cette méthode, on définit un ensemble de 
directions conjuguées par combinaison linéaire entre le gradient au point courant et la direction de recherche choisie 
à l'itération précédente. Cela guarantit la convergence du critère en au plus n itérations, n désignant la taille du 
problème traité. 

Dans le cas où le critère à minimiser n'est pas quadratique, nous utilisons une des méthodes suivantes : 

- l'algorithme du gradient conjugué non linéaire : cette méthode découle de la méthode du gradient conjugué linéaire 
et consiste donc à définir une direction de recherche par combinaison linéaire entre le gradient au point courant 
et la direction de recherche choisie à l'itération précédente. 

- l'algorithme L-BFGS : cette méthode peut être vue comme une généralisation de la méthode du gradient conjugué 
non linéaire. Pour un entier m choisi par l'utilisateur, elle consiste à effectuer une combinaison linéaire entre le 
gradient au point courant et les directions de recherche choisies lors des m itérations précédentes. 

Dans les deux cas, nous écartons d'emblée la méthode de plus forte pente qui s'avère généralement moins efficace. 

2.3.2 Le choix d'un pas de progression 

Pour les différents algorithmes envisagés, il est nécessaire de définir un pas de progression a à chaque itération après 
la définition de la direction de recherche d. Pour cela, on considère la fonction $ dont les variations sont celles du critère 
le long de la direction considérée : ^(a) — C(x + ad). A priori, le pas retenu doit correspondre à un minimum local : on 
doit rechercher a tel que <i>'(a) — 0. 
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Pour l'algorithme du gradient conjugué linéaire (cas de la minimisation d'un critère quadratique), le minimiseur du 
critère est unique et s'obtient de manière analytique. Lorsque le critère n'est pas quadratique, le minimiseur s'obtient 
de manière itérative. Dans ce cas, il est préférable de retenir un pas vérifiant les conditions dites de Wolfe [ ]. L'intérêt 
est double : d'une part, cela permet de s'approcher d'un minimum de <i> pour un nombre d'évaluations du critère et du 
gradient raisonnable. D'autre part, on s'assure d'obtenir un pas de progression suffisamment proche d'un minimum pour 
que les méthodes de type gradient proposées donnent une direction de descente et que l'algorithme converge globalement. 
Cette procédure de minimisation est résumée dans l'algorithme 1. 



Algorithme 1 Algorithme itératif - Procédure générale de minimisation 
Initialisation 
répéter 

Calcul du critère et du gradient en Xk 
Définition d'une direction de recherche dk 
répéter 

Choix d'un pas de progression a 

Calcul du critère et du gradient en % + ad 
jusqu'à Vérification des conditions de Wolfe 
Xk+i ^ Xk + akdk 
jusqu'à Convergence 



On distingue deux conditions de Wolfe : 
La première condition de Wolfe ou condition d'Armijo : 

Cix + ad)<C{x) + ciaVC{xfd (2.22) 

La seconde condition de Wolfe : 

VC(x + ad)^d > C2Q!VC(x)^d (2.23) 

ou la seconde condition de Wolfe forte : 

I VC(x + ad)^d| < |c2«VC(x)^d| (2.24) 

Les coefficients Ci et C2 intervenant dans les inégalités précédentes doivent être choisis tels que : < Ci < C2 < 1. 

Pour déterminer un pas de progression, nous utiliserons l'algorithme proposé par Moré et Thuente [ ] qui a l'avantage 
de déterminer un pas satisfaisant les conditions fortes de Wolfe pour un nombre limité de calculs du critère et du gradient. 

L'ensemble des méthodes présentées ici ont l'avantage d'être bien adaptées aux problèmes de grande taille puisque 
l'on définit une nouvelle direction de descente pour un coiit de calcul faible (on la calcule par simple combinaison linéaire 
des gradients aux itérations précédentes et au point courant, on ne passe pas par le calcul du Hessien du critère). C'est 
la raison pour laquelle certaines méthodes de type Newton-Kantorovitch ne seront pas utilisées ici : elles nécessitent 
l'inversion d'une matrice (inversion du Hessien du critère ou d'une forme approchée) et seraient alors très coûteuses en 
calcul. De plus, un petit nombre de vecteurs suffit à définir la nouvelle direction de recherche à chaque itération, ce qui 
réduit l'espace mémoire nécessaire pour stocker les variables utilisées. 

2.4 Utilisation d'autres variables 

Les équations établies précédemment font intervenir les contrastes Xp 6t Xs- Cependant, parmi les différentes méthodes 
que nous aborderons, certaines montrerons un problème de sensibilité du critère vis-à-vis des variations de Xp ^t Xs autour 
de valeurs élevées, ce qui a tendance à ralentir la convergence des algorithmes. Nous serons alors amenés à utiliser d'autres 
variables notées (7p et Us afin d'améliorer la sensibilité du critère. Elles sont choisies de sorte que de faibles variations de 
Œp et CTg induisent de fortes variations de Xp et Xs pour des valeurs de contraste élevées. 

Nous proposons plusieurs changements de variable. On les liste dans le tableau suivant en donnant leurs relations par 
rapport à Xp et Xs ainsi que les valeurs caractéristiques de la terre et du béton qui leurs sont associées. On donne également 
l'allure des fonctions associées aux différents changements de variables proposés sur la Figure 2.1. 
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Variables proposées 


Expressions de Xp Xs 

t!ll iUllULlUll Lie (^p t;L C g 


Valeurs caractéristiques 


Terre 


Béton 


(7p ^ Vp 
O's — l^s 


'2 2 

Xp — ^ ^p.o 

2 ? 
Xs — CTj. 


O^p.ToiTO — 'JUU 
Cs, Terre = 150 


<7p, Béton — 'iUUU 
Cs, Béton = 2200 


ap = \/vp 


Xp = (1/c^p)^ - Wp,o 

Xs = - <0 


CTp, Terre = 3, 3.10 
CTs,Terre = 6,7.10"^ 


CTp, Béton = 2, 5.10 
CTs, Béton 4,5.10"^ 


Up — In Vp 
(Ts = In Vs 


Xp = cxp 2ap - -u^ 
Xs = exp 2crs - Vs.o 


CTp, Terre = 5,7 
CTs, Terre = 5, 


Cp, Béton = 8,3 
"■s, Béton = 7, 7 




Figure 2.1 - Allure des fonctions associées aux différents changements de variable (dans cet exemple, Vp^ — Vg.o = 2) 



2.5 Données utilisées pour tester les différentes méthodes 

Afin de comparer les performances des méthodes abordées, nous avons effectué des tests sur un jeu de données syn- 
thétiques générées à l'aide de l'algorithme de résolution du problème direct. Nous avons travaillé sur un milieu de taille 
réduite puis sur un milieu de taille intermédiaire. Nous ne travaillerons pas sur des milieux de taille réelle (dimensions 
plus grandes et résolution plus fine) car la place mémoire requise serait trop importante. 

Nous présentons ci-dessous les deux milieux utilisés. Sur chaque figure, on représente : 

- en bleu clair la partie du milieu D qui reste invariante (le contraste reste nul, les caractéristiques restent celles du 
milieu de référence) ; 

- en rouge la zone d'étude (zone dans laquelle les caractéristiques évoluent), l'élément en béton est représenté en rouge 
plus foncé ; 

- en bleu foncé la zone PML ; 

- en jaune la position des capteurs (on mesure la composante verticale de la vitesse) ; 

- en vert les différentes positions de la source. 

Pour les différentes méthodes d'inversion proposées, seules les caractéristiques des pixels appartenant à la zone d'étude 
évoluent. Par conséquent, nous présenterons les résultats obtenus en n'affichant que le contenu de la zone d'étude. 

2.5.1 Milieu de petite taille 

Le milieu de petite taille entouré de la zone PML est représenté sur la Figure 2.2. 
On donne ci-dessous ses caractéristiques : 

- la taille du milieu est de 1 m de profondeur et 2 m de largeur ; 

- le signal source est un ricker centré à 200 Hz ; 

- on retient 15 fréquences réparties de façon uniforme entre 46,7 Hz et 700 Hz ; 

- la source est positionnée à 0,2 mètre en profondeur et est placée successivement à 0,5 m, 1 m et 1,5 m sur l'axe 
horizontal ; 
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10 20 30 40 50 60 

Figure 2.2 - Milieu utilisé pour les premiers tests avec les capteurs (jaune), les positions successives de la source (vert), 
la zone PML (bleu foncé), la zone d'étude contenant le bloc de béton (rouge) 



- les capteurs sont positionnés à 0,2 mètre en profondeur et sont espacés de 5 cm sur toute la longueur du milieu ; 

- le milieu recherché comprend un bloc de béton (20 cm de hauteur, 25 cm de largeur) entouré de terre (pas d'air) ; 

- la résolution en x et en y est de 0,05 m ; 

- il n'y a pas d'atténuation des ondes ; 

- la taille de la zone PML est de 50 cm. 

2.5.2 Milieu de taille intermédiaire 

On représente le milieu de taille intermédiaire entouré de la zone PML sur la Figure 2.3. 




50 100 150 200 260 300 



Figure 2.3 - Milieu de plus grande taille avec les capteurs (jaune), les positions successives de la source (vert), la zone 
PML (bleu foncé), la zone d'étude contenant le bloc de béton (rouge) 

On donne ci-dessous ses caractéristiques : 

- la taille du milieu est de 1 m de profondeur et 14 m de largeur ; 

- le signal source est un Ricker centré à 200 Hz ; 

- on retient 15 fréquences réparties de façon uniforme entre 46,7 Hz et 700 Hz ; 

- la source est positionnée à mètre en profondeur et est placée de à 14 m par pas de 0,7 m sur l'axe horizontal; 

- les capteurs sont positionnés à mètre en profondeur et sont espacés de 5 cm sur toute la longueur du milieu ; 

- le milieu recherché comprend est une superposition de dalles de béton (hauteur totale : 0,5 m, largeur totale : 1,15 
m) entourées de terre (pas d'air) ; 

- la résolution en x et en y est de 0,05 m ; 

- il n'y a pas d'atténuation des ondes ; 

- la taille de la zone PML est de 50 cm. 
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Chapitre 3 



Méthodes d'inversion fondées sur une 
formulation bilinéaire 

Toutes les méthodes basées sur une formulation bilinéaire suivent une démarche commune pour déterminer la solution 
du problème d'inversion. Nous allons détailler, dans une première partie, les spécificités de ce type d'inversion, puis nous 
présenterons chaque méthode ainsi que les résultats obtenus. 

3.1 Idée générale de l'inversion basée sur la formulation bilinéaire 

Le modèle direct permet, à partir des équations de propagation et d'une fonction source connue, de calculer les 
composantes du champ de vitesse en tout point d'un milieu souterrain dont on connaît les caractéristiques physiques. Dans 
le problème d'inversion, les caractéristiques physiques du milieu sont inconnues, de même que les composantes du champ 
de vitesse à l'intérieur du domaine. Seules sont connues les composantes verticales du champ de vitesse que l'on mesure en 
surface du milieu considéré. Dans le cas des formulations bilinéaires, on estime non seulement les variables d'intérêt liées 
aux caractéristiques du milieu, mais aussi celles liées au champ de vitesse, en tentant de compenser l'augmentation du 
nombre de paramètres à estimer par une plus grande simplicité de la formulation algébrique du problème. Ceci se traduit 
par des procédures d'estimation plus faciles à mettre en œuvre et plus efficaces algorithmiquement . De plus, nous avons 
le choix entre travailler avec des variables représentant les grandeurs physiques, ou utiliser le contraste de ces grandeurs 
par rapport à un milieu de référence. Ce choix permet de définir un certain nombre de variantes de la méthode d'inversion 
bilinéaire et a une influence directe sur les caractéristiques numériques des méthodes d'estimation. Dans ce qui suit, nous 
donnons une formulation générale du problème, et nous tentons de mettre en évidence les éléments qui ont l'impact le 
plus significatif sur les caractéristiques de la méthode d'inversion. 

3.1.1 Variables à estimer 

Comme indiqué précédemment, la formulation bilinéaire met en jeu deux ensembles de variables a; et ^ que l'on définit 
comme suit : 

Les variables d'intérêt Xp et Xs qui ne dépendent que de la position et qui paramétrisent les caractéristiques du milieu 
que sont les vitesses des ondes P et S 

Les variables auxiliaires z^^k qui dépendent de la fréquence et de la position de la source et qui paramétrisent le champ 
de vitesse V^^^k- 

Comme nous l'avons dit précédemment, les variables Xp, Xg et z^j^k peuvent être choisies de plusieurs manières, ce qui influe 
sur les caractéristiques des méthodes d'estimation correspondantes. De plus, il est possible d'introduire des changements 
de variables scalaires sur les quantités Xp, Xs et z^j.k- L'effet de ces choix est décrit au paragraphe 3.1.5. 

3.1.2 Système d'équations 

Quelle que soit la nature exacte des Xp, Xg et z^^ ^j la relation entre le milieu inconnu et les mesures î/^^^fc est modélisée 
à l'aide d'un système de deux équations qui prennent la forme générale suivante : 
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Une équation d'observation qui relie le vecteur des mesures yuj,k aux variables auxiliaires à travers la fonction TZ{z^^k) ■ 

Uuj.k = 'EiTZ{zuj,k) + buj,k 

où El est la matrice d'échantillonnage aux capteurs définie dans la section 2.1.3, b^j^k désigne un bruit qui représente 
les erreurs de mesure et de modélisation, et TZ{zui,k) est une fonction des variables auxiliaires dont l'expression sera 
détaillée ultérieurement. 

Une équation de couplage qui relie les variables auxiliaires aux variables d'intérêt et qui peut être considérée comme 
une contrainte sur le champ total de vitesse. Cette équation prend la forme générique suivante : 

■H(2„,fc) ='K{xp,Xs)'R{z^^k) 

oii Hi^Zi^^k) est une fonction de Zi^^k dont la forme sera précisée ultérieurement. 

La principale propriété du système d'équations ci-dessus est son caractère bilinéaire, ce qui signifie que chacune des 
équations est linéaire (ou plus précisément affine) par rapport aux variables Xp, Xg, et par rapport à z. En particulier, de 
par la linéarité de l'équation d'observation par rapport à z, on a nécessairement : 



où Md et Ud désignent respectivement une matrice et un vecteur constants. 
De même, l'équation de couplage est linéaire par rapport à 2, ce qui implique : 



(3.1) 



(3.2) 



avec Me une matrice constante et Uc un vecteur constant. Par ailleurs, la linéarité de l'équation de couplage par rapport 
à X permet d'affirmer que les éléments de la matrice K sont des combinaisons linéaires des Xp et Xg. De plus, on a vu 
dans la section 1.4 que la matrice d'impédance peut se décomposer en une somme de trois matrices indépendantes. De la 
même façon, la matrice K peut se décomposer en trois matrices K^', K* et K" (cette dernière pouvant être nulle suivant 
la méthode employée). et K'* sont respectivement des combinaisons linéaires des Xp pour la première et des Xs pour la 
deuxième. On a alors les deux égalités suivantes : 



Kî'7^(z^,fc) 
K^7^(z^,fe) 



^u,k^P 



,k^s 



(3.3) 



où j, et sont des matrices qui dépendent linéairement de z. Ces égalités permettent d'exprimer explicitement la 

linéarité des équations par rapport à chacun des jeux de variables. 

En utilisant la décomposition des matrices K^* et K*^ donnée dans la section 1.4, on établit les expressions explicites 
des matrices A^ et Af^ ^ : 



vP - 



■Diag{aî}H^ 
Diag{a^}Hy 

■Diag{a?}Hy 
Diag{aî}H'' 

Diag{a^}H^' 
Diag{af}Hy 

-Diag{aî}H'' 
-Diagja^'IHy 



Diag{[Diag{af}G- Diag{a|}Gy] 7^(^„,fc)} 
Diag{[Diag{a|}Gy Diag{af}G-] 7^(z^,fc)} 
-Diag{a|}Gy] 7^(z^,fe)} 



(3.4) 



Diag{[Diag{af}G^ 
Diag{[Diag{a|}G 



Diag{a|}Gyl 7^(z^,fe)} 



Dans les équations ci-dessus, les matrices G et H dépendent linéairement de z. Il est remarquable de constater que la 
structure des matrices A^ ^ et A^ est indépendante du choix précis des variables Xp, Xs et z. 



3.1.3 Critère 
Forme générale 

L'objectif de l'inversion est de déterminer les valeurs des différentes variables qui sont solution des équations de couplage 
et d'observation. Pour cela, nous allons construire un critère à partir d'une pénalisation quadratique de l'erreur sur chacune 



D. Vautrin, M. Voorons, J. Hier, Y. Goussard 



24 



Chapitre 3. Méthodes d'inversion fondées sur une formulation bihnéairc 



de ces deux équations. C'est en cherchant les valeurs des variables qui minimisent ce critère, et donc la somme des erreurs 
sur chacune des équations, que nous trouverons la solution au problème d'inversion. 

Minimiser l'erreur sur l'équation d'observation permet de s'assurer que la solution sera compatible avec les données 
mesurées par les capteurs, tandis que l'erreur sur l'équation de couplage doit être vue comme une contrainte que les 
variables doivent satisfaire pour que le modèle physique décrivant le problème soit respecté en tout point du domaine. 
Cette contrainte est relâchée, c'est-à-dire qu'on tolère que l'équation de couplage ne soit pas strictement vérifiée. Pour 
contrôler l'importance du respect de la contrainte et régir le compromis entre l'erreur qui touche l'équation de couplage 
et celle qui affecte l'équation d'observation, nous modifions le critère en introduisant une pondération de l'erreur sur 
l'équation de couplage par un hyperparamètre 7c. 

Pour faire face au caractère mal posé du problème, il est possible d'introduire des connaissances a priori sur la nature et 
le comportement spatial des différentes variables à inverser. Ceci peut être fait par le biais d'une fonction de régularisation 
qui est composée de plusieurs termes, chacun traduisant mathématiquement une information connue sur une des variables. 
Chaque terme est pondéré par un hyperparamètre qui rend compte de la confiance en l'information qu'il apporte. Le critère 
que nous utilisons dans le cadre des méthodes fondées sur une approche bilinéaire est donc constitué de 3 termes : 

OÙ : 

- Cd correspond à l'erreur de mesure ; son expression se déduit de l'équation d'observation 

Cd{xp,Xs,z^^k) = huiM-'EiTZiz^^kW 

- Ce correspond à l'erreur sur l'équation de couplage ; ce terme peut s'écrire de trois façons en faisant apparaître 
explicitement chaque ensemble de variables grâce aux jeux d'équations (3.1), (3.2) et (3.3) : 

Cc{xp,Xs,z^,k) = \\H{z^m) ~f<-{xp,Xs)TZ{z^MW 

= \\n{z^^k) ~ Al.xp - A^^feX, - K-7^(z„,fe)f 

= Il (Me - K{Xp,Xs)Md)Zu,.k + Me - K{Xp,Xs)ud\\^ 

- 7c est l'hyperparamètre qui quantifie le relâchement de la contrainte. Il peut dépendre des variables ou être constant. 

- (j) est la fonction de régularisation qui permet d'introduire de l'information a priori sur les variables Xp, Xs et Zi^^k- 

Régularisation employée 

La fonction de régularisation que nous avons employée pour toutes les formulations bilinéaires présentées est composée 
de trois termes : un rappel aux valeurs caractéristiques de la terre (xj et xj) sur les variables Xp et x^, un autre de rappel 
à zéro sur les variables z^j^k et un terme de pénalisation quadratique des différences premières de Xp et x^. 

Le terme de rappel à la valeur de la terre de la variable d'intérêt (qui devient un rappel à zéro dans le cas oîi l'on 
utilise des contrastes) correspond à la traduction d'une connaissance a priori sur la distribution spatiale de cette variable. 
Nous faisons l'hypothèse que la taille de l'objet difîractant enfoui est petite par rapport à celle de la zone d'étude, et donc 
que la plus grande partie des pixels correspondent à de la terre. Dans le cas de l'utilisation d'un contraste sur la variable 
d'intérêt, cela revient à dire que l'on considère que la plus grande partie de la zone d'étude est identique au milieu de 
référence et donc de contraste nul. 

Le terme de rappel à zéro sur les variables auxiliaires correspond à une connaissance a priori peu précise sur le champ 
de vitesse. Il pénalise les valeurs élevées du champ de vitesse, et donc permet d'éviter l'apparition de vitesses arbitrairement 
grandes ; il permet en outre d'améliorer le conditionnement de la méthode d'inversion. 

Le terme de pénalisation quadratique des différences premières a pour effet de pénaliser les grandes variations spatiales 
des variables d'intérêt et favorise ainsi l'apparition de zones homogènes. Toutefois, ce terme présente certaines limitations 
dans notre cas, puisque les valeurs des caractéristiques physiques recherchées peuvent varier fortement d'un objet à l'autre, 
et que pénaliser quadratiquement la variation spatiale de ces variables ne favorise pas la reconstruction d'une image aux 
transitions franches entre objets, introduisant ainsi une incertitude sur le positionnement des frontières. Cependant, les 
changements de variables proposés dans la section 2.4 permettent de réduire l'écart entre les valeurs des vitesses des ondes 
P et S des différents objets, et donc minimisent l'importance de cette affirmation. Pour mieux prendre en compte cet 
aspect du problème, nous avons envisagé de remplacer ce terme par une pénalisation de type L2L1, mais le temps nous a 
manqué pour mettre en œuvre et tester cette modification. 
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3.1.4 Gradients 

Toutes les méthodes d'optimisation envisagées pour résoudre le problème inverse sont décrites dans la section 2.3 et 
utilisent le gradient du critère. Dans le cas des méthodes fondées sur l'approche bilinéaire, la minimisation du critère (3.5) 
se fait suivant une direction de descente qui dépend du gradient par rapport aux variables d'intérêt et aux variables 
auxihaires. Cela nécessite donc d'étabhr l'expression du gradient du critère pour chaque jeu de variables. 

Tous calculs faits, le gradient relativement à z^^ ^. est donné par l'expression suivante : 

g,^ , {oj, k) = V,^_,C(2^,fe) = 2((M);EÎEiM<i + 7,(M, - Vi{xp, a:,)Mrf)^(M, - K{xp, x,)md))z^^k (3.6) 

- {m\¥.\{y^^k + Eiîtd) - 7,(M, - K(a;p, a;,)Mrf)^(M, - K(a;p, a;,)^^))) + V ,^ J{z^^k) 

Les gradients par rapport à Xp et Xg sont calculés en utilisant l'expression Ce faisant explicitement intervenir Xp et Xg. 
Ces deux gradients sont donnés par les équations suivantes : 

Q., = y.SlcCc{xp) + <l,{xp)) = 2^,{Y,Y.KkK,kXp - K\kin^^,k) - A^,,x,)) + V^Xxp) (3.7) 

k 

G.^ = y.AlcCcixs) + Hxs)) = 27,(^^ Af,^,A:,,,.x, - A^^,(H(z„,,) - Al .xp)) + V,>(x,) (3.8) 

k u 

L'expression des gradients de 4> par rapport à chacune des variables est détaillée dans la section 3.1.6. 

3.1.5 Choix des variables à inverser 

Les formulations bilinéaires du problème ne diffèrent les unes des autres que par les variables sur lesquelles on travaille. 
On a le choix de travailler directement sur les variables physiques, à savoir le carré des vitesses des ondes P et S ou les 
composantes du champ de vitesse, ou de travailler sur les contrastes et les sources de contrastes définis dans la section 2.1.2. 
Le choix des variables à inverser implique des différences techniques et numériques qu'il est important de souligner pour 
comprendre pleinement les avantages et inconvénients de chaque méthode. 

Implications techniques 

Dans le cas oii les Xp et Xs correspondent aux contrastes des variables d'intérêt par rapport à un milieu de référence, 
il est alors possible de restreindre l'estimation de cette variable à une zone d'étude comme nous l'avons signalé dans la 
section 2.1.5. En utilisant une connaissance a priori sur la position et la taille de l'objet enfoui, nous pouvons considérer 
que les caractéristiques physiques du sous-sol ne diffèrent de l'arrière-plan de référence que dans une petite partie du miheu 
D. Les contrastes ne prennent donc des valeurs non nulles que dans une zone restreinte de l'image, et ils ne doivent plus 
être calculés que dans la zone d'étude ce qui diminue d'autant le nombre d'inconnues à déterminer. La réduction de la zone 
de travail permet aussi de s'affranchir de la dépendance fréquentielle des matrices d'impédance puisqu'on fait l'hypothèse 
que les contrastes sont nuls dans la zone des PML. Il ne faut donc plus stocker en mémoire et mettre à jour qu'une seule 
matrice de contraste Xp et pour chaque caractéristique physique du sol, ce qui représente un gain intéressant en place 
mémoire et en temps de calcul. 

Par contre, si l'on choisit d'inverser le carré des vitesses des ondes P et S, il est alors nécessaire de faire l'estimation 
de ces variables en tout point du milieu D, incluant les PML. Cela nous oblige à mettre à jour et conserver en mémoire 
une matrice d'impédance par fréquence. 

A titre d'exemple, dans le cas du milieu de petite taille présenté à la section 2.5.1, on passe de 2604 inconnues à 150 
si on n'estime les contrastes que dans la zone d'étude suggérée. La matrice d'impédance est creuse et ne possède que 18 
lignes non nulles, comme nous pouvons le voir à la figure (1.3). Nous travaillons en double précision donc une matrice 
d'impédance occupe environ 6 Mb. L'utilisation des contrastes permet de ne stocker que deux matrices de 6 Mb contre 
une matrice 6 Mb par fréquence dans le cas où l'on utilise le carré des vitesses des ondes P et S (les résultats présentés 
dans ce document ont été obtenus en utilisant 15 fréquences). 

Le choix des variables auxiliaires à inverser est aussi un point important de la méthode d'inversion. Les variables 
auxiliaires z^^^]^ introduites dans la formulation bilinéaire dépendent de la fréquence et de la position de la source. Le 
nombre total de composantes de z^^j^k à déterminer est donc élevé, ce qui va demander un temps de calcul important. 
L'utilisation des sources de contrastes, définies dans la section 2.1.2, permet également d'en restreindre l'estimation à une 
zone d'étude, d'oii un gain important en nombre total d'inconnues et en volume de calcul. Par contre, si l'on choisit de 
travailler directement avec les composantes du champ des vitesses, on est obligé de les estimer en tout point du domaine 
D. Si on reprend le cas du domaine de petite taille, on passe de 300 x Nj x inconnues à déterminer si on utilise les 
sources de contraste sur une zone restreinte, à 5208 x Nf x TV^ dans le cas oîi on ne réduit pas la zone d'étude. 
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Un autre avantage de travailler sur une petite zone, outre les gains en temps de calcul et en place mémoire, est que 
cela a tendance à réduire la sous-détermination du problème. Effectivement, le nombre d'inconnues diminue tandis que le 
nombre de mesures ne change pas. Cependant, pour que l'utilisation de la zone d'étude soit justifiée, il faut disposer d'un 
milieu de référence qui soit assez proche de la réalité, ce qui reste jusqu'à présent un problème ouvert. 

Implications numériques 

L'utilisation de tel ou tel jeu de variable modifie les équations de couplage et d'observation, et définit donc une 
nouvelle méthode d'inversion. Le critère diffère suivant la méthode, et l'expression du gradient par rapport à l'une ou 
l'autre des variables fait intervenir des formes algébriques différentes plus ou moins simples et rapides à calculer. Le 
choix des variables à estimer a donc un impact direct sur la complexité calculatoire de l'inversion et le comportement 
numérique des méthodes d'optimisation. En particulier, les estimateurs associés à chacune des méthodes d'inversion ont 
des conditionnements différents, ce qui influe fortement sur la vitesse de convergence et la qualité du résultat. 

3.1.6 Détail de la mise en œuvre des algorithmes d'inversion 

Nous avons vu que le choix des variables à inverser a des conséquences importantes sur les caractéristiques numériques 
d'un algorithme d'inversion basé sur une formulation bilinéaire du problème. Les différents choix de couples de variables à 
inverser permettent de créer autant de variantes de la méthode d'inversion bilinéaire, dont la mise en œuvre nécessite un 
certain nombre de choix techniques tels que la forme de la fonction de régularisation utilisée, les stratégies d'optimisation 
possibles, les changements de variables envisagés, la prise en compte progressive des fréquences dans l'inversion, la déter- 
mination des hyperparamètres, les données utilisées pour les tests, les critères de convergence employés et les problèmes de 
conditionnement de la matrice d'impédance du problème direct. Nous allons discuter dans cette section des divers aspects 
numériques et techniques soulevés par ces choix de mise en œuvre. 

Détail de la régularisation employée 

Le premier choix à faire concerne la forme de la fonction de régularisation. De celle-ci va dépendre l'expression des 
estimateurs de chaque variable, ce qui va influer sur le type de méthode d'optimisation à employer. Nous avons vu dans 
la section 3.L3 que la fonction de régularisation était composée de trois termes. Cette fonction, qui sera utilisée dans 
toutes les méthodes d'inversion basées sur une approche bilinéaire abordées dans ce document, est donnée par l'expression 
suivante : 

<P{x,, xs,z) = 7;o(l|2;p - {^vhf + \\xs - {xshf) + 7^0 E E ll^-.fcll' + 7."i(ll(Di + ^2)x,r + ||(Di + T>2)xsf) (3.9) 

k UJ 

OÙ la notation (• )t fait référence aux caractéristiques de la terre, Di et D2 sont les matrices des différences premières 
dans les directions verticale et horizontale, et les 7^q, 7^q et 7^]^ sont les hyperparamètres qui pondèrent l'importance des 
différents a ■priori introduits par la régularisation (7^0 pour le rappel à la terre des variables d'intérêt, 7^q pour le rappel à 
zéro sur les variables auxiliaires et pour la pénalisation quadratique des différences premières des variables d'intérêt). 

La fonction de régularisation (3.9) est quadratique par rapport à chacune des variables. Le gradient d'une telle fonction 
par rapport à chacune des variables se calcule simplement. On obtient : 

V,^0(a:p) = 2-f^.^{xp ~ (2:^)7) + 27,^i(Di + D2)*(Di -|- Ds)^^ (3.10) 
V^Mx.) = 2^r^{xs - {xs)t) + 27."i(Di + D2)*(Di + Ds)^^ (3.11) 
V.^^,0(z„,fe) = 27^0 Y. E (3.12) 

k LO 

Stratégie d'optimisation employée 

Pour réaliser une inversion basée sur une formulation bilinéaire, nous pouvons choisir d'optimiser les deux ensembles de 
variables soit de façon alternée, soit simultanément. De plus, il est possible d'estimer conjointement, ou de façon alternée, 
les différentes variables de chaque jeu de variables. Plus précisément, l'estimation des deux variables d'intérêt Xp et Xg peut 
être fait alternativement ou conjointement, et il en va de même pour les variables auxiliaires z^j.k correspondant à chaque 
fréquence et chaque position de tir. 

Nous avons choisi d'inverser les deux variables d'intérêt conjointement. Les valeurs de Xp et Xg sont du même ordre de 
grandeur, et leur estimation respective partage un certain nombre de calculs, comme on peut le voir sur les équations (3.7) 
et (3.8). Il est donc avantageux de réaliser l'estimation de ces deux variables simultanément. Il suffit pour cela de concaténer 
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les systèmes linéaires de chacun des estimateurs. Cette stratégie a été appliquée à toutes les méthodes présentées dans ce 
chapitre. 

De même, il existe deux façons de voir l'estimation des variables auxiliaires z^j^k- Soit on effectue successivement 
la résolution d'autant de systèmes linéaires qu'il y a de positions de la source et de fréquences, soit on estime toutes les 
fréquences conjointement en un seul bloc. Des tests nous ont montré que les deux approches donnent des résultats similaires. 
Nous avons donc choisi d'utiliser l'algorithme d'inversion alternée qui est détaillé dans le pseudo-code de l'algorithme (2). 



L'estimation alternée des deux jeux de variables repose sur la construction itérative et alternée de séquences (fj, ^p)'-"-' 
et {(4j,/c)*'"''; w € [1 . . . Nf] et k € [1 . . . N^]}- Les variables sont estimées alternativement en utilisant une méthode d'op- 
timisation tronquée, chaque inconnue est estimée avec seulement quelques itérations de l'algorithme d'optimisation. Une 



Algorithme 2 Algorithme d'inversion alternée 



Entrées: Initialisation des valeurs Xp^\ a;i°^ et zj^l. 
n = l 
répéter 

pour LU = LOi, . . . ,LûNf faire {Boucle sur le nombre de fréquences Nf} 

pour k = 1, . . . , Nk faire {Boucle sur le nombre de positions de la source N^} 

Déterminer (fcj.fe)'-"-' G argmin C{zi^^k) à l'aide de (f^;,*:)^"""'^-' 

fin pour 
fin pour 

Déterminer £p"^ € argmin C{xp) à l'aide de £p" 



Déterminer fi" G argmin C{xs) à l'aide de fi 

Mise à jour des matrices d'impédances 
n = n + l 
jusqu'à Convergence 



résolution complète serait coiiteuse en temps de calcul et n'est pas forcément utile dans le cadre d'une optimisation par 
blocs. 

Comme nous l'avons vu dans la section 2.3, divers algorithmes d'optimisation sont disponibles pour résoudre un système 
linéaire. Nous avons constaté que les résultats obtenus à l'aide de LBFGSB et du gradient conjugué sont similaires, mais 
que le gradient conjugué est généralement plus rapide. Nous avons donc choisi d'utiliser ce dernier pour résoudre les 
systèmes linéaires intervenant dans l'estimation de chaque variable. 

La méthode d'optimisation alternée peut s'avérer inefficace si le minimum de la fonction de coût se trouve dans une 
vallée très étroite. Dans ce cas, la recherche du minimum global en suivant alternativement la direction donnée par le 
gradient associé à chaque jeu de variables se fait selon une trajectoire en zigzag, et un nombre très important de petits pas 
est requis pour l'atteindre [( >]. Pour remédier à ce problème, il est possible de réaliser l'optimisation simultanée de toutes 
les variables, c'est-à-dire de construire itérativement une suite d'estimées (fp, fs, {4j,fc; ^ & [l . . . Nf] et k £ [1 . . . iVfc]})^"^ 
De cette façon, on n'effectue la descente du gradient que suivant une seule direction ; on limite ainsi les problèmes de 
zigzag liés à la méthode alternée. L'optimisation simultanée des deux ensembles de variables est réalisée en concaténant 
les différents systèmes algébriques. 

Le pseudo code de la méthode est donné dans l'algorithme (3). 

Lorsque l'on réalise l'estimation simultanée des Xp, Xg et Zi^^k, le pas de mise à jour de la méthode d'optimisation ne 
s'obtient plus en résolvant un système linéaire [ ]. Comme indiqué dans la section 2.3.2, nous utilisons l'algorithme de 
Moré Thuente pour calculer le pas de descente et la méthode d'optimisation choisie est le LBFGSB qui semble converger 
plus rapidement que la méthode du gradient conjugué non linéaire. L'inconvénient de cette méthode d'inversion simultanée 
est qu'il faut s'assurer que toutes les variables inversées simultanément ont la même échelle pour que la convergence ne 
soit pas trop lente. 



Influence des changements de variables 
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Algorithme 3 Algorithme d'inversion conjointe 
Entrées: Initialisation des valeurs x^, a;° et 2° 

n = 1 

répéter 

Déterminer (fp, fs, 4;,^)'^"^ e argmin C{xp, Xs, z^^k) à l'aide de {xp, Xg, iui^kY"^'^^ 

Mise à jour des matrices d'impédances 
n = n + \ 
jusqu'à Convergence 



Les changements de variables proposés dans la section 2.4 ont pour particularité de modifier significativement l'intervalle 
des valeurs que peuvent prendre Xp et Xg, tout en n'entraînant qu'une petite modification du calcul du gradient. On tente 
ainsi de modifier le conditionnement du problème sans augmenter significativement le volume de calcul. 

Chaque changement de variable de x en cr étant de nature scalaire, la modification apportée au calcul du gradient 
consiste simplement à multiplier ce dernier par une matrice diagonale : 

V<,C(a) = Diag(V<,a;)V,C(a;) 



Introduction progressive des fréquences 

Les méthodes présentées précédemment réalisent l'inversion en tenant compte de toutes les fréquences simultanément, 
comme dans [ ]. Cependant, les travaux de Pratt et al. [9, lU, 11] ont montré qu'il pouvait être intéressant d'utiliser une 
stratégie d'incorporation progressive des fréquences pour effectuer l'inversion. Les basses fréquences ont un comportement 
plus linéaire vis-à-vis des perturbations du modèle que les hautes fréquences [10]. Il est donc suggéré par Pratt [11] 
d'effectuer l'inversion en partant d'un certain nombre de basses fréquences et d'introduire, au fur et à mesure de l'inversion, 
de nouvelles fréquences plus élevées pour ajouter l'information qu'elles contiennent. Cette technique a aussi l'avantage 
d'accélérer les premières itérations en réduisant considérablement le nombre de variables à estimer. 

L'algorithme (4) décrit la méthode utilisée pour introduire les fréquences progressivement. Les fréquences sont intro- 
duites par paquets de deux, des plus basses vers les plus hautes. 



Algorithme 4 Algorithme d'inversion avec introduction progressive des fréquences 
Entrées: Initialisation de la méthode d'inversion 

Entrées: Regroupement des Nf fréquences en Nq groupes Gi de 2 fréquences, i E [l,No] 
G = 

pour i = 1, . . . , Nq faire {Boucle sur les groupes de fréquences, des plus basses fréquences vers les plus hautes} 
n = 1 

Construction du groupe de fréquences utilisé pour l'inversion G = G U Gi 
répéter 

Estimation des contrastes x avec une itération de la méthode d'inversion et les fréquences du groupe G 
Calcule du déplacement relatif de la solution var = ""'•'"(x^-Xn-i) 

^ norm{Xn — l) 

n = n + l 
jusqu'à n > rimax ou var < seuil 
fin pour 



Détermination des hyperparamètres 

Les hyperparamètres sont tous fixés empiriquement. Ils pondèrent l'importance de la régularisation et de la contrainte 
imposée par l'équation de couplage vis-à-vis de l'adéquation aux données. 

Données utilisées 

Tous les résultats présentés dans ce chapitre ont été obtenus en utilisant des données synthétiques. Les mesures sont 
issues du modèle direct développé par EdF auxquelles on a ajouté un bruit blanc gaussien tel que le rapport signal sur bruit 
soit de 30 dB. La géométrie et les caractéristiques du milieu utilisées pour les générer sont décrites dans la section 2.5.1. 
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Nous nous sommes servi de deux jeux de données pour initialiser les méthodes lors des tests. Le premier jeu est 
obtenu en utilisant les caractéristiques physiques de la terre, et le deuxième correspond à la solution exacte du problème. 
L'initialisation à la solution permet d'apprécier la dégradation introduite par la méthode d'inversion et sert de référence. 

Critère de convergence employé 

Le critère d'arrêt que nous avons choisi d'employer est un seuil maximal sur la valeur de la norme du gradient qui 
est fixé empiriquement. Toutefois, nous utilisons également le déplacement relatif de l'estimée pour stopper l'algorithme. 
Cette métrique mesure la variation de l'estimateur considéré entre les itérations n et n — 1, et permet de voir si l'estimée 
que l'on calcule à chaque itération stagne à cause de problèmes numériques ou continue à évoluer. Par exemple, dans le 
cas du déplacement relatif de la variable Xp, il est défini comme : 



Ikp-'ll 

Problèmes de conditionnement de la matrice d'impédance du problème direct 

La matrice d'impédance, définie dans la section L3, est construite en utilisant les équations physiques (1.7) décrivant 
le problème de propagation des ondes dans un milieu élastique. Ces équations lient les inconnues du problème entre elles, 
les valeurs des vitesses Vp et Vg aux composantes du champ de vitesse Vx et Vy. Ce lien est algébriquement réalisé par 
l'introduction d'une matrice d'impédance. A partir des équations de la physique, nous avons formulé une méthode pour 
générer cette matrice en utilisant une décomposition de l'opérateur de discrétisation en plusieurs matrices assimilables à 
des filtres (voir partie L4). 

La construction de la matrice d'impédance du modèle direct, telle que mise en œuvre initialement par EdF et décrite 
dans [2], est basée sur un domaine (composé de la zone d'étude et des PML) ceinturé par une bande de 1 pixel de large 
de valeur 1 autour du domaine. Ceci se traduit, pour une image de taille M x N, par l'ajout de 2{M + N ~ 2) valeurs 
égales à —1 sur la diagonale de la matrice d'impédance. Cet ajout n'affecte pas le problème direct et donne des résultats 
similaires à ceux obtenus en utilisant la matrice d'impédance construite tel que nous le présentons dans la section L4. 
L'erreur quadratique moyenne relative entre les résultats du problème direct générés en utilisant la matrice d'impédance 
que nous proposons et celle fournie par EdF est de 10"^"* en moyenne sur toutes les fréquences. 

Toutefois, nous avons observé des différences significatives au niveau algébrique entre les deux matrices d'impédance. 
La présence de ces valeurs supplémentaires égales à —1 sur la diagonale de la matrice fournie par EdF est responsable d'une 
importante dégradation du conditionnement de la matrice d'impédance A^^. Une estimation du nombre de conditionnement 
des matrices calculées par le code fourni par EdF nous indique qu'il est de l'ordre de 10^^ (en fixant les variables d'intérêt 
à la solution). Cependant, en utilisant les matrices d'impédances construites en utilisant la méthode montrée dans la 
section L4, le nombre de conditionnement de A(^_p^s tombe à 10^. 

Une autre remarque, de moindre importance, peut être faite sur la structure des matrices d'impédance. Des collabo- 
rateurs au projet, travaillant chez EdF, nous ont suggéré d'effectuer la permutation des lignes et colonnes de ces matrices 
pour faire des gains de temps lors de certaines opérations. Notamment, si l'on permute les lignes et colonnes de ces matrices 
de telle sorte que l'on alterne les composantes et Vy de la vitesse dans le vecteur V^^k, on obtient des matrices d'im- 
pédance dont les valeurs non nulles sont concentrées autour de la diagonale, ce qui permet une accélération significative 
du temps de calcul d'un préconditionneur de type décomposition de Cholesky tronquée (10 fois plus rapide). Nous avons 
donc effectué le réarrangement des lignes et colonnes de ces matrices à chaque fois qu'un tel préconditionnement était 
employé. 

Les valeurs du nombre de conditionnement données dans ce chapitre sont toutes estimées à l'aide de la fonction condest 
de Matlab. 

3.2 Méthode CSI 
3.2.1 Description 

L'obtention des équations utilisées dans la méthode "Contraste Source Inversion" (CSI), [12], est détaillée dans la 
partie 2.1. La CSI est la méthode bilinéaire correspondant au choix de variables suivant : 

Les variables d'intérêt x sont les contrastes Xp 6t Xs du carré des vitesses Vp et Vg qui sont estimés sur une partie 
réduite du domaine D, et ont été définis dans la section 2.1.2 

Les variables auxiliaires z^ ^ sont les sources de contraste W^^k, définies dans la section 2.1.2, qui sont estimées sur la 
même partie réduite du domaine D 
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Nous rappelons ici les équations d'observation et de couplage, présentées dans la section 2.1, qui sont à la base de la 
méthode CSI : 

Équation d'observation y^^k = Ei V^ f. - Bf^ W^^k 
Équation de couplage Wu:,k = (Xp + X,)(E2 V^^^ - 

oîi les matrices BJ^ et Bf, sont définies dans la section 2.1.6 et E2 est la matrice d'échantillonnage sur la zone d'étude 
définie dans la section 2.1.5. 

En reprenant les notations de la section 3.1, on établit que : 

n{ W^^k) = W^^k avec =1 et «c = 

7^( W^^k) = fe - (A„,p.s)ô'E^ W^,k avec = -(A^,p,,)-iE* et - V^^^ (3-13) 

K(xp,Xs) = (Xp + X,)E2 

où {Ai^^p^s)q^ est l'inverse de la matrice d'impédance du modèle direct lié aux caractéristiques du domaine de référence 
qui a été définie dans la section 2.1. 

A partir de ce système d'équations, et en nous basant sur ce qui a été dit dans la section 3.1.3, nous définissons la 
fonction de coiit de la méthode CSI comme suit : 

k LJ 

+ 7c E E II ^"■'^ - (^P + ^^)(E2 V°^, - W^^k)f + (3.14) 

k 

À partir de l'expression du critère CSI (3.14), et en substituant les variables à inverser spécifiées par la méthode CSI 
et les notations données par (3.13), nous pouvons réécrire les formules du gradient par rapport aux variables d'intérêt et 
aux variables auxiliaires données par les équations (3.6), (3.7) et (3.8). L'expression du gradient du critère par rapport à 
chaque variable est donné ci-dessous : 

- Gradient par rapport à W^i^k ■ 

Qw^,, = V w^.,C{ W^^k) =2((B^tB^ + ^,(/ + (Xp + X,)B^)t(/ + (X^ + X,)B^) + 7.^) W^,^) 

+ B'^Jiv^.k - El 1/° ,) - 7c(/ + (Xp + X,)B^)t (Xp + X,)E2 V^^) 

- Gradient par rapport à Xp '■ 

Qx. - ^xAXp) =2(E T.^ilc^tkK,k + 7?o + 7?i(Di + D2)'(Di + J^,))xp 

k i>J 

-7cA^^,(W^.,fe-A^^,X.))) 

oîi et Af^ sont telles que : 

K.kXp = Xp[l/^% - (A„,p,,)o-iE* 
et Al^kXs - X,[I/0, - (A^,p,,)o-iE* 

Pour obtenir l'expression du gradient par rapport à Xs, il suffit d'intervertir les indices p et s. 

3.2.2 Spécificités de la méthodes CSI 
Aspects numériques 

On remarque que la matrice normale associée au calcul des W^^^k fait intervenir les matrices B"^ et B'^ qui sont des 
versions échantillonnées de (Aj^^p^^)^^. (A(^_p_s)o est une matrice creuse et son inversion est coiiteuse en temps de calcul et 
en place mémoire, surtout dans le cas d'un problème de dimension réelle. Nous avons décidé de ne pas inverser ces matrices, 
mais plutôt de résoudre les systèmes linéaires qui interviennent dans le calcul du gradient à l'aide de décompositions LU. 

La méthode CSI permet une estimation rapide des contrastes, d'autant plus qu'il est possible de calculer la matrice 
normale et d'utiliser un préconditionneur de Jacobi pour accélérer le calcul. 

Aspects techniques 
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La méthode CSI nous permet de restreindre l'estimation des variables Xp, Xs et Wuj,k à une petite zone d'étude tant 
que le milieu de référence est identique à l'arrière-plan de l'objet à détecter. L'estimation des sources de contraste s'avère 
être très lourde en temps de calcul. Pour chaque itération de la méthode d'inversion, il faut résoudre Nf x NgX {2 + Nac) 
systèmes linéaires, où Nf est le nombre de fréquences, Ng le nombre de positions de tir et Ncc le nombre d'itérations du 
gradient conjugué. 

A chaque itération, il faut mettre à jour les matrices j,, ^, Xp et X^. 
Le facteur de pondération CSI 

Les articles de Abubakar et Van Den Berg [12, L3, 14] suggèrent d'utihser une valeur de l'hyperparamètre telle que 
les deux composantes de la fonction de coût soient égales lorsque les sources de contraste sont nulles. On obtient ainsi 
l'expression suivante pour la valeur de l'hyperparamètre 7c : 

"^ = """ = E.EJI(Xp + x.)K.",.|P ^'-''^ 

Il s'est avéré que l'utilisation de ce facteur de pondération n'est pas adaptée à notre problème. Nous avons donc opté pour 
une estimation empirique des valeurs des hyperparamètres, ce qui nécessite des réglages supplémentaires fastidieux. 



3.2.3 Résultats 

Méthode d'optimisation alternée 

Les résultats que nous allons présenter ont été obtenus en utilisant les données décrites dans la section 3.1.6, et, pour 
des raisons de temps de calcul, nous nous sommes limité au milieu de petite taille. 

La figure (3.1) montre l'évolution du critère pour les différentes changement de variable proposés et pour deux ini- 
tialisations différentes. On constate que, quel que soit le changement de variable, le critère décroît très lentement. Même 
après un grand nombre d'itérations, aucune des méthodes ne semble pouvoir rejoindre la courbe du critère obtenue en 
initialisant l'algorithme avec la solution. 



Évolution du critère de ia méthiode CSI 



_ Init. Solution - x^ 

Init. Terre -x 

Init. Terre - 1/v 

Init. Terre - 1/\r 
p: 

Init. Terre - In v 

P.. 

init. Terre - v 

p.s 




1 1.5 
Itérations 



2.5 
xlO* 



Figure 3.1 - Évolution du critère de la méthode CSI pour les différents changements de variables et deux initialisations 



Les courbes de la figure (3.2) représentent l'évolution du déplacement relatif de l'estimée du contraste Xp &u cours des 
itérations. Une courbe semblable est obtenue dans le cas des Xs- On remarque que, assez rapidement, le déplacement du 
contraste d'une itération à l'autre est très petit, et qu'il ne cesse de décroître. L'évolution des variables se fait donc très 
lentement et de plus en plus lentement, ce qui explique pourquoi le critère décroit très lentement et semble stagner après 
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un certain nombre d'itérations. Cependant, le déplacement de la solution reste toujours suffisamment significatif pour ne 
pas fanchir le seuil fixé comme critère d'arrêt de l'agorithme tel que décrit dans la section 3.1.6. Il ne nous a donc pas été 
possible de faire converger cette méthode en un temps raisonnable. 



Évolution de la différence relative de % pour ia méttiode CSi 




Figure 3.2 - Évolution du déplacement relatif de l'estimée des contrastes pour les différents changements de variables et 
deux initialisations 



Les résultats obtenus pour la méthode CSI avec une optimisation alternée des deux jeux de variables, et pour les 
changements de variables indiqués à la section 2.4, sont présentés dans la figure (3.3). Les cartes montrées dans cette 
figure ont été obtenues après 17600 itérations (soit environ 3 jours de calcul). On constate que, même après un grand 
nombre d'itérations, aucune des méthodes testées ne permet de reconstruire l'amplitude et la forme de l'objet recherché. 
On peut noter toutefois, que tous les changements de variable testés ont un effet différent sur la convergence de la méthode. 
Notamment, le changement de variable en ln(î;p^s) donne les meilleurs résultats en termes de décroissance du critère et de 
différence relative, et les cartes obtenues semblent visuellement supérieures. 

En raison de la lenteur de l'évolution de la solution, il nous est impossible de dire si la méthode converge véritablement. 
La convergence des méthodes de type descente de gradient dépend fortement du nombre de conditionnement du hessien 
[15]. Or, on voit sur la figure (3.5) que le nombre de conditionnement du hessien des estimateurs CSI augmente au 
fil des itérations. Cette détérioration du conditionnement pourrait tout à fait expliquer les problèmes de stagnation du 
déplacement relatif des estimées et du critère mis en évidence dans les figures (3.1) et (3.2). Il semble donc légitime de 
penser que la lenteur de la convergence de la méthode CSI est due à un mauvais conditionnement des matrices normales 
associées aux estimateurs. 



Méthode d'optimisation conjointe 

Pour tenter de remédier aux problèmes de convergence de la méthode CSI alternée, nous avons mis en oeuvre une 
version d'estimation conjointe de la CSI, voir [.]. 

La mise en œuvre de cette méthode fait apparaître un problème de différence d'échelle entre les deux jeux d'inconnues, 
c'est-à-dire entre les contrastes et les sources de contrastes ; l'amplitude des premiers est de l'ordre de 10^, tandis que pour 
les deuxièmes elle ne dépasse pas 10^. Cela nuit au calcul du pas de descente utilisé par la méthode d'optimisation, qui sera 
un compromis entre les pas optimaux de chacune des deux variables. On résout traditionnellement ce type de problème 
en utilisant un préconditionneur de Jacobi ou un préconditionneur diagonal empirique basé sur une connaissance a priori 
de l'ordre de grandeur de chaque jeu de variables. Toutefois, les matrices normales des estimateurs de la méthode CSI ne 
sont pas à diagonale dominante, et il est reconnu que les préconditionneurs diagonaux ne sont pas efficaces dans ce cas 
[16]. Nous avons d'ailleurs effectué des tests qui nous ont convaincu que les préconditionneurs de Jacobi et empiriques ne 
sont pas adaptés à notre problème, car ils n'améliorent pas les résultats. 
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Les cartes de contraste calculées avec la méthode CSI conjointe sont présentées dans les figures (3.6) et (3.7). Elles 
ont été obtenues en pondérant le terme d'adéquation aux données du critère par la norme des mesures. Cette pondération 
permet de réduire l'importance du problème d'échelle entre les deux jeux de variables. Sans elle, le pas de progression de 
l'algorithme d'optimisation est trop petit, et l'évolution de la solution est excessivement lente. 

Les résultats montrent que cette méthode ne permet pas de régler les problèmes de lenteur de la convergence, elle 
semble même converger beaucoup plus lentement que la CSI classique. Les contrastes montrés dans le figure (3.6) sont 
qualitativement et quantitativement moins bons que ceux de la CSI après un temps de calcul équivalent. Cette méthode 
atteint plus rapidement le palier de stagnation de la solution, et ce, pour tous les changements de variables testés. On en 
déduit donc que, dans ce cas-ci, les meilleurs résultats sont obtenus sans cfi'cctuer de changement de variable. 



3.3 Méthode du gradient modifié 

Pour pallier les problèmes de conditionnement qui surviennent dans les méthodes de type CSI, nous avons exploré 
d'autres approches basées sur des formulations bilinéaires. Nous avons commencé par étudier la méthode du gradient 
modifié (CM). Elle est antérieure à la CSI et a été proposée la première fois par Kleinman et van den Berg [17]. 

3.3.1 Description 

Dans cette méthode, la variable d'intérêt est le contraste des caractéristiques physiques, et la variable auxiliaire repré- 
sente directement les composantes du champ de vitesse. Plus précisément, on a : 

Les variables d'intérêt x sont les contrastes Xp 6t Xs des carrés des vitesses Vp et Vs qui sont estimés sur une partie 
réduite du domaine D, et ont été définis dans la section 2.1.2 

Les variables auxiliaires z^^^k sont les composantes du champ de vitesse V^^j^k qui doivent être estimées sur tout le 
domaine D 

L'adéquation du modèle aux mesures et la relation bilinéaire entre les deux ensembles de variables, les contrastes Xp 
et Xs et les variables auxiliaires sont données dans le système d'équations suivant : 

Équation d'observation î/c^.fe = Ei 

Équation de couplage {A^^p^s)q{ V^ ^, - K;,^) = (Xp -|- X,,) V^^k 
En reprenant les notations de la section 3.1, on établit que : 

-RiV^.k) = (A^,p,,)o(KS,fc - V^^k) avec = -(A^,p,«)o et ît, = (A„,p,,)o ^ 

niV^.k) = V^^k avec M<j = 1 et ît^ = (3.16) 

K(xp,Xs) = Xp + Xs 

En se référant à la section 3.1.3 et en utilisant les équations de couplage et d'observation propres à la méthode CM, 
nous établissons que le critère prend la forme : 

k (jJ 

+ 7cEEll(Ac.,p.s)o(K.%- K.,fc)-(Xp + X,)V;,^fc||2 + </, (3.18) 

k w 

A partir de l'expression du critère, nous pouvons réécrire les formules des gradients donnés dans les équations (3.6), 
(3.7) et (3.8) en y substituant les équations (3.16), et en tenant compte du choix des variables fait par cette méthode. 
- Le gradient du critère par rapport à V^^^k est ; 

5v^,, = V v^,,C{ K;,fc) =2((7,(Xp + X, + (A„.p,,)o)^(Xp + X, + (A^,p,,)o) + EÎEi + 7^^^) V^^k 

- {{E\y^,k + 7c((Xp + X, + (A^,p,,)o)^(A<.,p,,)o F° ,)))) 
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- Le gradient du critère par rapport à Xp est : 

k u: 

avec 

L'expression du gradient du critère de la méthode GM par rapport à Xs est obtenu en interchangeant les indices p 
et s. L'estimation des Xp et Xs est faite conjointement. 

3.3.2 Spécificités de la méthode GM 

Le calcul des Xp et Xs est similaire au cas de la CSI ; seule l'expression de % change. 

Le calcul du gradient par rapport au champ de vitesse ne fait intervenir l'inverse d'aucune matrice et nécessite seulement 
l'évaluation de sommes de matrices et de produits matrice- vecteur. Le choix du champ de vitesse comme variable auxiliaire 
nous impose l'inversion d'un grand nombre de variables puisqu'on ne peut se restreindre à une zone d'étude pour la 
détermination de cette variable. Cet inconvénient est toutefois compensé par la rapidité et la simplicité du calcul du 
gradient par rapport aux variables auxiliaires. 

Lorsque l'on teste isolément l'estimateur des composantes du champ de vitesse de la méthode GM, en initialisant la 
méthode avec les contrastes solution, on constate que la convergence est extrêmement lente. Cet estimateur fait intervenir 
la matrice normale suivante : D^^ = A^^ ^ s-^^^,p,s + E|Ei où A(^_p g = (Xp + + (A^^ ^ s)o) est la matrice d'impédance 
du problème direct. Nous avons vu dans la section 3.L6 que cette matrice est très mal conditionnée. Le nombre de 
conditionnement associé à la matrice normale de l'estimateur du champ de vitesse de la méthode GM, en initialisant les 
contrastes à la solution, est de 10^^ lorsqu'on la génère avec le code fourni par EdF. Si on utilise la méthode suggérée dans 
la section L4 pour la construire, le conditionnement tombe à lO^'^. Même si cela reste élevé, il est maintenant possible de 
résoudre le système intervenant dans l'estimateur du champ de vitesse à l'aide d'une décomposition LU, ou d'utiliser des 
méthodes de préconditionnement basées sur une factorisation. 

Le choix d'un préconditionneur pour l'estimation du champ de vitesse est limité par le fait que les matrices normales 
associées ne sont pas à diagonale dominante. Les préconditionneurs diagonaux, comme celui de Jacobi, ou, quand ils 
existent, de type factorisation simple sont reconnus pour ne pas être efficaces pour ce type de problème [l i.]. Nous avons 
testé les préconditionneurs de Jacobi, LU incomplète et Cholesky incomplète, et nous avons remarqué qu'ils n'apportaient 
aucune amélioration de la convergence et pouvaient même parfois nuire. 

Il est cependant possible d'utiliser un préconditionneur issu de la décomposition de Cholesky incomplète décalée [18] 
qui est basé sur la factorisation de la matrice normale augmentée D^^ + al, où a est un paramètre assurant la stabilité de la 
factorisation qui doit être fixé manuellement et que l'on doit garder le plus petit possible. L'emploi de ce préconditionneur 
permet d'accélérer significativement la convergence des premières itérations de la méthode GM, mais l'évolution du critère 
finit par stagner. 

Il existe une autre formulation du gradient modifié qui ne diffère de celle qu'on vient de présenter que dans son équation 
de couplage qui prend la forme suivante ■ V° — V^^k = {^ui,p,s)Q^i^p + '^s)Vui,k- Cette variante de la méthode GM 
s'obtient facilement à partir de la méthode présentée en multipliant les deux membres de l'équation de couplage par 
{■^uj,p,s)^^ ■ La présence de la matrice (Atj^p^s)Q ^ dans l'expression des matrices normales des deux estimateurs ne permet 
pas de faire le calcul explicite du Hessien, et la complexité calculatoire de cette méthode est la même que celle de la CSI. 
Le calcul du gradient par rapport au champ de vitesses nécessite la résolution d'autant de systèmes linéaires que dans le 
cas de la méthode CSI, mais il doit être fait pour tous les points du domaine contrairement à la CSI. Pour ces raisons, 
nous avons choisi de ne pas tester cette méthode. 

3.3.3 Résultats 

Pour des raisons de temps de calcul, nous ne présentons les résulats de la méthode GM que pour le milieu de petit taille 
décrit dans la section 3.1.6. Les résultats ont été obtenus en utilisant une inversion directe par décomposition LU pour 
l'estimation du champ de vitesse et un gradient conjugué non linéaire pour l'estimation des contrastes. L'utilisation de 
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la décomposition LU entraîne un gain de temps important par rapport à l'emploi d'un gradient conjugué préconditionné 
avec un préconditionneur tel que suggéré dans la section 3.3.2. 

Les résultats de la méthode GM après 6000 itérations sont montrés dans les figures (3.8) et (3.9). On observe que 
le critère de la méthode GM décroit très lentement, et que l'utilisation des changements de variables décrits dans la 
section 2.4 influe peu sur cette situation. 

Cependant, le changement de variable en ln(wp_s) donne la meilleure décroissance du critère, mais les cartes de contraste 
obtenues sans changement de variable semblent plus intéressantes, car c'est dans ces conditions que l'erreur quadratique 
moyenne par rapport aux contrastes de la solution est la plus faible. 

La figure (3.10) montre le conditionnement de la matrice normale de l'estimateur du champ de vitesse. On remarque 
que celui-ci se dégrade fortement après seulement quelques itérations. Le système est donc mal conditionné, ce qui pourrait 
expliquer la lenteur de convergence de la méthode GM. 

La méthode GM a un comportement analogue à celui observé pour la méthode CSI en ce qui concerne la décroissance 
du critère, l'évolution du déplacement de la solution et le conditionnement des estimateurs. La méthode GM a donc le 
même problème de convergence que la méthode CSI ; elle s'avère néanmoins plus intéressante que cette dernière puisque le 
coût de calcul d'une itération est deux fois plus faible et que l'erreur quadratique moyenne atteinte après 17600 itérations 
est plus petite (avec 90.32% pour Xp et 82.27% pour Xs contre respectivement 93.56% et 90.43% dans le cas de la méthode 
CSI). 

3.4 Méthode sans contraste 

La méthode bilinéaire que nous présentons dans cette section ne fait pas appel à un milieu de référence et repose sur 
l'estimation directe des caractéristiques du milieu. Pour cette raison nous l'avons appelée méthode sans contraste. Nous 
avons étudié cette méthode car elle permet de faire intervenir des formes algébriques simples dans le calcul des gradients. 

3.4.1 Description 

Cette méthode travaille directement sur les grandeurs physiques. Plus précisément : 

Les variables d'intérêt Xp et Xs sont les carrés des vitesses des ondes P et S, Vp et V^, qui sont estimés sur tout le 
domaine D 

Les variables auxiliaires z^^k sont les composantes du champ de vitesse Kj.fc qui doivent être estimées sur tout le 
domaine D 

L'équation de couplage utilisée dans ce modèle est exactement celle utilisée par le modèle direct, et l'équation d'obser- 
vation correspond à un simple rééchantillonnage du champ des vitesses. 

Équation d'observation y^j.k — 'EiVi^^k 
Équation de couplage F„ & = A.^,p,s K^.fc 

où -Fcj.fe est un vecteur formé par les composantes de la fonction source, tel que défini dans la section 1.3. 
En reprenant les notations de la section 3.1, on établit que : 

■H( Vcj^k) = F^,k avec Me = et Mc = F^^k 

n{ V^^k) = K.,/c avec =1 et = (3.19) 

^(Xp: Xs) — ^u),p,s 

En se référant à la section 3.1.3 et à partir des équations de couplage et d'observation propres à la méthode sans 
contraste, nous établissons le critère comme : 

k oj k oj 

A partir de l'expression du critère, nous pouvons réécrire les formules des gradients donnés dans les équations (3.6), 
(3.7) et (3.8) en y substituant les équations (3.19) et en tenant compte du choix des variables fait par cette méthode. 
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- Le gradient du critère par rapport à V^i^k est : 

- Le gradient du critère par rapport à Xp est : 

- ^xAXp) -2(7c E E((^-fe^-fe + ^ro + 7?i(Di + D2)^(Di + D^)) V, 

k 

avec 

On obtient une expression du gradient semblable pour Vs, il suffit d'échanger les indices p et s. L'estimation des 
et Vg est faite conjointement. 

3.5 Spécificités de la méthode sans contraste 

Aspects numériques 

La méthode sans contraste est basée sur deux estimateurs très simples, et le calcul des gradients ne fait intervenir aucune 
résolution de système linéaire. Dans les deux cas, le Hessien peut être calculé, ce qui permet d'utiliser des méthodes de 
préconditionnement simples. Cependant, toutes les inconnues doivent être estimées sur tout le domaine, y compris les 
PML. L'estimateur du champ de vitesse est exactement le même que pour la méthode CM. Cette méthode souffre donc 
des mêmes problèmes de conditionnement. 

Aspects techniques 

La méthode d'inversion bilinéaire la plus simple ne fait pas intervenir de milieu de référence. Elle consiste à travailler 
directement sur les variables d'intérêt et à utiliser les composantes du champs de vitesse V^j^k comme ensemble de variables 
auxiliaires. L'avantage de cette méthode est que l'on s'affranchit de la difficulté de trouver un domaine de référence valable. 
Par contre, il faut estimer les deux jeux de variables sur tout le domaine, car elles sont a priori connues en aucun point 
du domaine, y compris dans les PML. Les matrices Ap et As dépendent donc de la fréquence, et il faut stocker et mettre 
à jour à chaque itération deux matrices A^^^p et A^^^g pour chaque fréquence. 

3.5.1 Résultats 

L'estimateur des Vg et Vp, lorsque testés séparemment en fixant les V^^^k à la solution, se comporte bien et converge 
vers la solution en une centaine d'itérations du gradient conjugué, si l'on utilise un préconditionneur de Jacobi. Sans le 
préconditionneur de Jacobi, il faut compter autour de 100 fois plus d'itérations pour que cet estimateur converge vers la 
solution. Le critère de convergence atteint est la décroissance de la norme de la solution qui est inférieure à un seuil, et 
l'erreur quadratique moyenne obtenue après convergence est de l'ordre de 10~^. 

L'estimateur Vui,k converge très lentement. La matrice AIj p ^A^^p^s + e|Ei, qui intervient dans le calcul de cet 
estimateur, est mal conditionnée. Son nombre de conditionnement est de l'ordre de 10^^. Un préconditionnement simple 
basé sur la diagonale du Hessien ne permet pas de résoudre le problème de mauvais conditionnement. Cependant, le 
problème de conditionnement peut être traité en utilisant la même approche que celle exposée dans la partie 3.3. Dans ce 
cas, l'estimateur des champs de vitesses converge rapidement vers la solution. 

Les résultats présentés dans la figure (3.11) correspondent aux itérations 21000 pour l'initialisation à la terre et 
8200 pour l'initialisation à la solution. Le calcul des variables d'intérêt se fait sur tout le domaine, donc les cartes des 
caractéristiques de la terre sont données pour tout le domaine sauf les PML. 

Même si les deux estimateurs se comportent bien lorsqu'ils sont testés séparément, les résultats obtenus par cette 
méthode sont moins intéressants que ceux obtenus avec la méthode CSI. Après 21000 itérations, l'amplitude maximale du 
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signal reconstruit n'est que de 2, 5 x 10^ dans le cas de Xp st 1.6 x 10^ pour Xs, alors qu'on atteignait 2 x 10® et 1, 5 x 10® 
dans le cas de la CSI. Il faut cependant signaler que le calcul d'une itération est très rapide, deux fois plus rapide que pour 
la CSI. On constate également, voir figure (3.12), que cette méthode souffre du même problème de lenteur de convergence 
que toutes les autres. 

3.6 La méthode CFSI 

La méthode CFSI [ i ] (contraste field source inversion) est une version augmentée de la CSI dans laquelle on relâche la 
contrainte sur la définition des sources de contrastes, W^^k = (Xp + Xs) Vi^,k, en incluant le terme d'erreur correspondant 
dans l'expression du critère. On introduit ainsi un jeu de variables auxiliaires supplémentaire : les composantes du champ 
de vitesses. L'estimation de ces composantes ne peut être retreinte à une zone d'étude, car, à l'exception du tour du 
domaine sur lequel elles doivent satisfaire les contraintes de bord de l'équation (1.2), elles ont une valeur non nulle en tout 
point du domaine. Du fait de l'introduction d'une nouvelle contrainte et d'un nouveau jeu de variables, cette méthode ne 
rentre pas tout à fait dans le cadre général défini dans la section 3.1, mais l'approche employée pour résoudre le problème 
et les choix techniques de la mise en œuvre sont analogues. 

3.6.1 Description 

Trois jeux de variables doivent être estimés : 

Les variables d'intérêt Xp et Xs sont les contrastes Xp et Xs des vitesses et qui sont estimés sur une partie réduite 
du domaine D, et ont été définis dans la section 2.1.2. 

Les variables auxiliaires z^j^k sont les sources de contraste qui sont estimées sur une zone restreinte du domaine D, et 
ont été définies dans la section 2.1.2 

Les variables auxiliaires supplémentaires sont les composantes du champ de vitesses V^i^k qui sont à estimer sur 
tout le domaine D 

La CFSI repose sur trois jeux d'équations qui se déduisent directement de la formulation CSI, l'adéquations aux 
données et deux contraintes : 

Équation d'observation y^^k = Ei V^ ,. - B^J W^^k 

Équation de couplage V^^k = V° f. - B^^ Wu;,k 

Contrainte sur les sources de contraste FFcj.fc = O^p + '^s)Vuj.k 

Cette méthode diffère des autres méthodes bilinéaires par l'introduction d'un quatrième terme dans l'expression du 
critère qui correspond à l'erreur quadratique sur la définition des sources de contraste. A partir de ces équations on forme 
donc le critère suivant : 

C = EE ll?'-.^- - [El V°,k - BtiW^,k)]f 

k cj 

+-fcJ2Y.\\V"^'^~^°.k+'Bzw^,kr 

k uj 

+ 7«. E E II ^-.fe - + ^-.'=11' + 

k uj 

Cette méthode entraîne l'introduction d'un nouvel hyperparamètre 7^1 servant à pondérer l'importance de l'erreur sur 
l'équation des sources de contraste. Comme pour les hyperparamètres des méthodes précédentes, celui-ci est estimé empi- 
riquement. 

La minimisation du critère de la CFSI par rapport à chaque variable nécessite de formuler l'expression des gradients 
par rapport à chacun des trois jeux d'inconnues : 
- Gradient du critère par rapport à Vi^^k '■ 

Sv^,, = V v^,,Ci K.,fe) = 2((7cl + 7^(Xp + X,)t (Xp + X,)) V^^k 

- (7.(Xp + X,)t W^^k + Ici fe - W^^k))) 
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- Gradient du critère par rapport à Wcj,k ■ 

- (7»(Xp + X,) - 7cB^n - - Bt}{y^^, - Ei 1^° ,))) 



- Gradient du critère par rapport à Xp '■ 

k Lû 

avec 

L'expression du gradient du critère de la méthode CFSI par rapport à Xs est obtenu en interchangeant les indices p 
et s. L'estimation des Xp et Xs est faite conjointement. 

3.6.2 Spécificités de la méthode CFSI 
Aspects numériques 

L'intérêt de la CFSI par rapport à la CSI est qu'elle permet d'avoir une structure algébrique du calcul des WL,fe 
différente. Cependant, cela nécessite l'introduction de nouvelles variables, ce qui ajoute un nombre important d'inconnues 
au problème. On se trouve ainsi devoir estimer les contrastes et sources de contraste sur la zone d'intérêt, et les composantes 
du champ de vitesse sur tout le domaine. 

Il est néanmoins important de souligner que l'estimateur V^^^k a pour avantage de ne faire intervenir l'inverse d'aucune 
matrice. Le hessien peut donc être calculé facilement, ce qui permet de résoudre le système à l'aide d'une décomposition 
LU, ou d'utiliser un algorithme de type gradient conjugué préconditionné. Le nouveau jeu de variables introduit par la 
méthode CFSI n'ajoute donc pas une charge de calcul trop importante puisque celui-ci peut s'estimer de façon relativement 
simple et rapide. 

La complexité du calcul de W^^^ est sensiblement la même que pour la CSI, mais fait intervenir des structures 
algébriques différentes. L'estimateur des contrastes de la CFSI est très semblable à celui de la méthode CSI, sauf pour ce 
qui est du calcul des matrices j, et ^ . 

Le principal inconvénient de la CFSI est qu'elle introduit un hyperparamètre supplémentaire qu'il est difficile de 
fixer empiriquement. Si on teste séparemment chaque estimateur, c'est-à-dire en initialisant avec la solution des deux 
autres variables, alors W^i.k, Xp et Xs convergent rapidement vers la solution. Cependant, le conditionnement des matrices 
normales de ces estimateurs dépend directement du choix des hyperparamètres. Par exemple, le calcul de V^,k fait intervenir 
la matrice normale 7Jl-t-7tu(Xp-|-Xs)^(Xp-|-Xs) dont le conditionnement est de l'ordre de 10^° si on utilise les contrastes 
solution pour construire les matrices X^ et Xp et que les hyperparamètres sont égaux à 1. Dans le même cas de figure, la 
matrice normale de M^j.fc est, quant-à elle, bien conditionnée. Les courbes de la figure (3.13) montrent le logarithme du 
conditionnement des matrices normales des estimateurs Vu,,k et W^j^k en fonction des valeurs des hyperparamètres 7c et 
7^. On peut voir que le conditionnement de ces deux matrices varie fortement suivant la valeur des hyperparamètres. 

D'une manière générale, il est préférable que 7^, soit inférieur à 7c. Le choix des ces paramètres résulte d'un compromis 
entre le conditionnement des deux matrices normales. Les deux ne peuvent être simultanément bien conditionnées, il faut 
donc choisir les hyperparamètres pour que le conditionnement de chacune des matrices ne soit pas trop mauvais et que 
la méthode ne diverge pas au bout d'un certain nombre d'itérations. Les tests empiriques ont montré que l'intervalle de 
valeurs qu'ils peuvent prendre est relativement étroit et qu'en dehors de cet intervalle la méthode diverge ou converge 
excessivement lentement. Si 7c est choisi dans l'intervalle [10~^, 1] et 7^ dans [10^^^, 10^^°], alors la méthode semble 
converger. 

L'utilisation d'un rappel à zéro comme fonction de régularisation permet, dans une certaine mesure, de réduire le 
problème de conditionnement. Cependant, le poids que l'on peut donner à ce terme de rappel à zéro est limité par la 
dégradation qu'il introduit sur la solution. L'amélioration du conditionnement ainsi obtenue est donc limitée. 
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3.6.3 Résultats 

Les résultats obtenus pour la méthode CFSI après 1380 itérations sont présentés aux figures (3.14) et (3.15). 

Les cartes de la figure (3.14) montrent qui existe un certain nombre de pixels possédant une forte amplitude négative 
qui rendent difficile l'appréciation des résultats. 

Une fois de plus, nous constatons que l'évolution de la solution est lente. De plus, le temps de calcul d'une itération 
est deux fois plus important que pour la CSI. 

3.6.4 Inversion en multifréquence 

Dans la section 3.1.6, nous avons décrit une stratégie d'incorporation progressive des fréquences pour réaliser l'inversion. 
Nous avons testé cette stratégie sur les méthodes CSI et CM sans utiliser de changement de variable. Les résultats sont 
présentés aux figures (3.16) et (3.17). Comme les méthodes testées convergent très lentement, il était exclu d'attendre la 
convergence de la méthode pour un nombre de fréquences donné avant de procéder à l'introduction de nouvelles fréquences 
dans l'algorithme d'inversion. Nous avons choisi d'incorporer une fréquence à chaque fois que la différence relative des 
contrastes serait inférieure à un certain seuil, tout en assurant un minimum de 2000 itérations. 

Les courbes de la figure (3.18) montrent le déplacement relatif des contrastes. On remarque que l'introduction d'un 
nouveau groupe de fréquences entraîne un déplacement de la solution plus important pour les basses fréquences que 
pour les hautes fréquences. Ceci se traduit par une accélération de la convergence initiale ; cependant, la méthode semble 
atteindre un palier à partir duquel son comportement est similaire à celui observé lors de l'inversion simultanée de toutes 
les fréquences. 
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(a) Initialisation à la solution, variable Xp,j 




(b) Initialisation aux caractéristiques de la terre, variable Xp,^ 





(c) Changement de variable en l/vp^s, initialisation aux caractéristiques de la terre 




(d) Changement de variable en l/fp s, initialisation aux caractéristiques de la terre 



3? 




(e) Changement de variable en In(t)p^s), initialisation aux caractéristiques de la terre 

Figure 3.3 - Résultats des contrastes Xs (à droite) et Xp (à gauche) obtenus avec la méthode CSI pour différents 
changements de variables et deux initialisations 
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(a) Changement de variable en fp,s, initialisation aux caractéristiques de la terre 



Figure 3.4 ~ Résultats des contrastes Xs (à droite) et Xp (à gauche) obtenus avec la méthode CSI pour différents 
changements de variables et deux initialisations (suite) 




(a) Évolution du conditionnement de l'estimateur des(b) Évolution du conditionnement de l'estimateur des 
variables auxiliaires pour plusieurs fréquences contrastes 



Figure 3.5 - Évolution du conditionnement des estimateurs de la méthode CSI 
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(b) Initialisation aux caractéristiques de la terre, variable Xp.s 




(c) Changement de variable en l/fp,s, initialisation aux caractéristiques de la terre 




(d) Changement de variable en log(vp_s), initialisation aux caractéristiques de la terre 

Figure 3.6 - Résultats de la méthode CSI conjointe 
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Évolution du oritère de la méthode CSI conjointe 



3 4 
Itérations 



. Init. Solution -x^ 
Init. Terre -y 

Init. Terre - 1/v 

p,s 

Init. Terre - In v 
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Figure 3.7 - Évolution du critère de la méthode CSI conjointe 
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(a) Initialisation à la solution, variable Xp,j 
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(b) Initialisation à la terre, variable Xp,. 




(c) Changement de variable en l/vp^s, initialisation aux caractéristiques de la terre 




X 10 

I 



I 




X 10 
10 



(d) Changement de variable en l/ii^ ^, initialisation aux caractéristiques de la terre 




I 



I 



X 10 

3 

2.5 
2 

1.5 



I 



(e) Changement de variable en ln(ï)p_s), initialisation aux caractéristiques de la terre 

Figure 3.8 - Résultats obtenus avec la méthode GM pour différents changements de variables et deux initialisations 
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Évolution du critère de la méthode GM 



Init. Solution - % 
Init. Terre -X^^ 
Init. Terre - 1/v 
Init. Terre - 1A 
_ Init. Terre - In 



p,s 



3000 
Itérations 



Figure 3.9 - Évolution du critère de la méthode GM pour différents changements de variables et deux initialisations 
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(a) Évolution du conditionnement de l'estimateur de(b) Évolution du conditionnement de l'estimateur des 
la variable auxiliaire pour plusieurs fréquences contrastes 



Figure 3.10 - Évolution du conditionnement des estimateurs de la méthode GM en initialisant à la solution et à la terre 



D. Vautrin, M. Voorons, J. Idier, Y. Goussard 



46 



Chapitre 3. Méthodes d'inversion fondées sur une formulation bihnéairc 





(a) Initialisation à la solution, variable Xp,s 





(b) Initialisation à la terre, variable Xp,s 

Figure 3.11 - Résultats obtenus avec la méthode sans contraste 



Evolution du critère de la méttiode basique 



_ Init. Solution -x^ 
_ Init. Terre -x 



1 1.5 
Itérations 



Figure 3.12 - Évolution du critère de la méthode sans contraste 
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Conditionnement de i'estimateur de V 




(a) Conditionnement de V^^k 



Conditionnement de i'estimateur de W 




(b) Conditionnement de W^.fc 

Figure 3.13 - Logarithme du conditionnement des matrices normales en fonction des hyperparamètres 
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(a) Initialisation à la solution, variable Xp,s 



r 



(b) Initialisation aux caractéristiques de la terre, variable Xp,^ 

Figure 3.14 - Résultats de la méthode CFSI 



Evolution du critère de la méthode CFSI 
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_ Init. Terre -% 
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Figure 3.15 - Évolution du critère de la méthode CFSI 
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(a) Initialisation avec la solution 
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(b) Initialisation à la terre 



Évolution du critère 
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2000 4000 6000 8000 10000 12000 14000 16000 

Itérations 

(c) Évolution du critère 

Figure 3.16 - Résultats obtenus par la méthode CSI en introduisant progressivement les fréquences 
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(a) Initialisation avec la solution 
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(b) Initialisation à la terre 
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Evolution du critère 
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(c) Évolution du critère 

Figure 3.17 - Résultats obtenus par la méthode GM en introduisant progressivement les fréquences 
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Evolution du contraste en P 



E 10 




Évolution du contraste en S 
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(a) Évolution du déplacement relatif du contraste Xp 




2000 4000 6000 8000 10000 12000 14000 16000 
Itérations 

(b) Évolution du déplacement relatif du contraste Xs 



Figure 3.18 - Évolution du déplacement relatif des contrastes obtenus par la méthode GM en introduisant progressivement 
les fréquences 
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Méthodes inverses fondées sur la 
formulation primale du problème direct 

Dans le cas où l'on ne considère que le terme d'adéquation aux données, le critère associé à la formulation primale 
s'écrit de la façon suivante : 

C{Xp,Xs) = -^^^hi^^k - 9uj,k{Xp,XsW (4-1) 

w k 

OÙ : 

- w et fc désignent respectivement la fréquence et la position de la source considérées. 
~ yu:,k désigne le jeu de données mesurées ; 

~ 5w,fc(XpiXs) correspond à la formulation primale du problème direct et désigne la fonction de construction des 
données synthétiques pour une distribution de contrastes donnée ; 

9c.AXp,Xs) = (V^^°fc)capt-Bd(I+(Xp + X,)B^)-i(Xp + X,)X.% (4.2) 
(la construction de cette formulation est détaillée dans le Chapitre 2, page 15) 

4.1 Les procédures itératives envisagées 

Pour minimiser le critère, nous envisageons d'utiliser une procédure itérative consistant à chaque itération à déterminer 
une direction de recherche dans l'espace de représentation puis à chercher un pas de progression efficace le long de cette 
direction de recherche (minimisation en une dimension). 

En ce qui concerne la définition d'une direction de recherche, nous avons choisi de tester les algorithmes du gradient 
conjugué non linéaire et L-BFGS tel que décrit dans la Partie 2.3 page 18 (le critère minimisé n'est pas quadratique). 

Nous avons également envisagé de tester une méthode d'optimisation pixel par pixel. Cette méthode consiste à optimiser 
une des grandeurs d'intérêt (xp ou Xs) en un pixel, les autres caractéristiques restant constantes. Cela revient à considérer 
les variations du critère parallèlement à l'un des axes de l'espace d'état ; on ne prend donc pas en compte les variations 
locales du critère pour choisir la direction de descente. Malgré le fait qu'un nombre important de balayages de l'image soit 
nécessaire, cette méthode s'est avérée efficace pour certains problèmes d'imagerie par ondes diffractées car l'optimisation 
des caractéristiques en un pixel est très rapide [10]. 

Concernant le choix d'un pas de progression, les valeurs des coefficients intervenant dans les conditions de Wolfe sont 
généralement telles que ci ~ 10"'' et C2 — 0, 1 avec l'algorithme du gradient conjugué non linéaire [ ]. Nous avons retenu 
ces valeurs pour les tests, les résultats seront présentés par la suite. 

4.2 Proposition d'une formulation primale optimisée 

Avec la formulation primale actuelle (Eq. 2.21, page 18), le calcul du critère (terme d'adéquation aux données) et de 
son gradient en un point quelconque de l'espace de représentation nécessite la résolution de deux systèmes linéaires pour 
chaque fréquence et chaque position de la source. Les matrices normales de ces sytèmes sont égales à I + (Xp + Xs)B^. 
Elles dépendent de la fréquence et des contrastes mais ne dépendent pas de la position de la source. Nous proposons donc 
de passer par la décomposition LU de ces matrices afin de diminuer le coût de calcul. 
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Cependant, la matrice est une matrice pleine. Les matrices normales des systèmes linéaires sont donc pleines, ce 
qui nécessite une quantité d'espace mémoire importante. On s'attend donc à ce que la résolution des systèmes linéaires 
(décomposition LU des matrices et résolution des systèmes linéaires triangulaires qui s'ensuivent) soit longue. De plus, les 
matrices qui interviennent dans l'expression du critère et du gradient sont également pleines. 

Nous proposons d'utiliser une formulation primale optimisée ne faisant intervenir que des matrices creuses. Nous 
allégerons ainsi le temps et l'espace mémoire requis pour effectuer les calculs. 

4.2.1 Modifications apportées à la formulation primale 



Nouvelle expression faisant intervenir une matrice normale creuse 



Nous commençons par réécrire l'expression de la matrice normale des systèmes linéaires en faisant intervenir l'inverse 



de 



(I + (Xp + X,)BS)-i = (B^)-i((B^)-i + Xp + X,)-i (4.3) 

La matrice (B^)^^ est une matrice creuse. En effet, on sait que B^ est une sous-matrice carrée de AJq (voir Partie 
2.1.6, page 17) et que A^^^o est creuse. A une permutation des lignes et des colonnes près, décomposons AJq et A(^^o en 



quatre sous-matrices : 



(A.,o)-' = 



B^ 



et A^^o = 



où Mi^uj, M2,uj, M'i,^ et A^4,tj sont creuses et Mi^uj est de la taille de B^. 
Appliquons ensuite le lemme d'inversion par bloc ^ . On obtient : 

(B^)-i - Mi,^ - M2,^MllM3,. 



(4.4) 



(4.5) 



Si l'on désigne par nzE,x et nzE,y les dimensions horizontales et verticales de la zone d'étude en nombre de pixels : 

- Mi^uj est une matrice creuse de taille 2{nzE,x + ^){nzE,y + 1) x 2{nzE,x + ^){nzE.y + 1) (chaque ligne contient au 
plus 18 éléments non nuls) ; 

- M^Ij est une matrice pleine ; 

- -^2,1^ est une matrice creuse. Elle compte 2(nzE,x + l)('^ZE,y + 1) lignes. Seules 4:{nzE,x + nzE,y) lignes de cette 
matrice contiennent des éléments non nuls. Il s'agit du nombre de composantes de -Fi^.fc qui s'expriment en fonction 
de composantes de V° n'appartenant pas à la zone d'étude (voir équation 2.1 page 15) ; 

- Ms,ùj est une matrice creuse. Elle compte 2(nzE,a: + ^)(nzE,y + 1) colonnes. Seules 'i{nzE,x + nzE,y) colonnes de 
cette matrice contiennent des éléments non nuls. Il s'agit du nombre de composantes de ^ appartenant à la zone 
d'étude et intervenant dans les expressions des composantes de F^^^k qui n'appartiennent pas à la zone d'étude. 

Si nzE,a; et nzE,y sont suffisamment élevés, on en déduit que la matrice Al2,w-A^4 i-A^s.w contient une majorité de 
coefficients nuls. La matrice (B^)^^ est donc creuse. 

En tenant compte de cette nouvelle expression, le calcul du critère et du gradient passe maintenant par la résolution 
de systèmes linéaires dont les matrices normales sont creuses. On procédera donc de la manière suivante : 

- lors d'une phase d'initialisation, on calculera les matrices (B^)^^ associées aux différentes fréquences en utilisant la 
décomposition de la matrice A^^^o en quatre sous-matrices ; 

- pour chaque point de l'espace de représentation considéré, on calculera la matrice (B^)^^+Xp+Xs puis on effectuera 
sa décomposition LU ; 

- pour résoudre les systèmes linéaires intervenant dans le calcul du critère et du gradient, on utilisera les facteurs 



et U 



(B<= )-i+Xp+Xs s-fin de se ramener à la résolution de systèmes linéaires triangulaires. 
Comme la matrice (B^)"^ -|- Xp -|- Xg ne dépend pas de la position de la source, cette procédure est plus économe que 
l'inversion de chaque système pris indépendamment. 



Modification supplémentaire 

Les matrices BjJ qui interviennent dans les expressions du critère et du gradient sont des matrices pleines. Ce sont 
également des sous-matrices de A~q (voir Partie 2.1.6, page 17). On propose de les remplacer par leur expression faisant 
intervenir Af 



'1 



- Ei(A„_p^s)o ^E*2 



(4.6) 



1. Lemme d'inversion par bloc : 



'A 


B 


-1 y 


C 


D 





A-1 + A-'^B(D - CA-'^B)-^CA-^ 
-(D- CA-i_B)-iCA-i 



-A-iB(D- CA-iS)-i 
(D - CA-^B)-^ 
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Ainsi, au lieu de multiplications par B^, nous effectuons des inversions de systèmes linéaires dont la matrice normale 
est creuse. Etant donné que la matrice A^,.o ne dépend pas du contraste, les mêmes matrices normales interviendront à 
plusieurs reprises. On procède donc de la manière suivante : 

- lors d'une phase d'initialisation, on effectue la décomposition LU des matrices A(^_o ; 

- pour chaque système linéaire associé, on utilise les facteurs ï^a^^o Ua„,o '^^ ramener à la résolution de 
systèmes linéaires triangulaires. 

Cela permet de gagner notamment en place mémoire et en temps de calcul lors de la phase d'initialisation (matrices 
longues à calculer). 

4.2.2 Nouvelle écriture de la formulation primale 

Désormais, nous écrirons la relation guj,k{XpjXs) Qui lie les contrastes aux données synthétiques de la façon suivante : 

gu:AXp,Xs) = «Jcapt - Ei(A„,p,,)o-^E|(BS)-iMji(Xp + X3)K.% (4.7) 
avec : M„ = (B^)"! + Xp + 

4.2.3 Mise en évidence du gain en temps de calcul et en espace mémoire 

L'utilisation de la formulation optimisée permet un gain à la fois en temps de calcul et en espace mémoire occupé. 
Nous le montrons en utilisant les expressions avant et après modifications sur le cas suivant : 

- le domaine considéré est constitué de 3700 pixels et la zone d'étude comprend 800 pixels ; 

- le nombre de fréquences retenues est égal à 15 et la source prend trois positions différentes. 

Mise en évidence du gain en temps de calcul 

Nous présentons dans le tableau suivant le temps de calcul associé aux étapes les plus coûteuses pour le calcul du 
critère et du gradient en un point. Les temps de calculs affichés correspondent à la prise en compte de l'ensemble des 
fréquences et des positions de la source. 



Avec la formulation initiale 


Avec la formulation optimisée 


Initialisation 

Calcul des Bj, et des 


1 heure 


Initialisation 

Calcul des (B^)"^ 
Décomposition LU des A^^^q 


35 secondes 

30 secondes 
5 secondes 


Calcul du critère et du gradient 

Calcul des I + (Xp -f Xs)Bj 
Décomposition LU des I + (Xp + Xs)B^ 
Calcul du critère 

(1 sysicmc linéaire par fréq. et par pos.) 
Calcul du gradient 

(1 système linéaire par fréq. et par pos.) 


> 1 minute 

10 secondes 
50 secondes 
2 secondes 

2 secondes 


Calcul du critère et du gradient 

Décomposition LU des (B^)"""" + Xp + Xg 
Calcul du critère 

(2 systèmes linéaires par fréq. et par pos.) 
Calcul du gradient 

(2 systèmes linéaires par fréq. et par pos.) 


6 secondes 

2 secondes 
2 secondes 

2 secondes 



Mise en évidence du gain en mémoire 

Nous mettons en évidence le gain en espace mémoire obtenu en utilisant les expressions avant et après modification 
sur le même domaine (nous n'affichons que les éléments occupant le plus d'espace mémoire) : 



Avec la formulation initiale 


> 120 Mo / frcq 


Avec la formulation optimisée 


- 21 Mo / frcq 




40 Mo / freq 




7 Mo / freq 




1 Mo / freq 




7 Mo / freq 


I + (Xp + Xs)B= 


40 Mo / frcq 




1.5 Mo / freq 


Ll+(Xp+Xs)Bj 


20 Mo / freq 


(B^)-i+Xp + X, 


1.5 Mo / freq 


Ul+(Xp+Xs)BC 


20 Mo / freq 


'^(Bjj-l+Xp+Xa 


2 Mo / freq 






^(B?,)-1+Xp+Xs 


2 Mo / freq 
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4.3 Calcul du gradient du critère 

Pour les différentes méthodes d'optimisation proposées, il est nécessaire de savoir calculer le gradient en un point 
quelconque de l'espace d'état. Nous résumons ici une démarche permettant d'aboutir à l'expression du gradient du terme 
d'adéquation aux données en un point quelconque. Pour le terme de régularisation, il n'y a pas de difficulté. 

Expression du gradient du critère pour une seule fréquence et une seule source 

Nous nous plaçons dans le cas où l'on ne considère qu'une seule fréquence et qu'une seule position de la source afin de 
simplifier le raisonnement et les écritures. On a donc : 

C{Xp,Xs)^\\\y~9{Xp.XsW avec g(xp, X.) = Kapt " ^iAq-'E* (B'^)-iM-i(Xp + X,)T/0 (4.8) 

où M = (B'=)"i +Xp +Xs. 

Les gradients selon Xp et Xs sont égaux à : 

dC BC 

= -nMxp,Xs)\y ~ g{xp,Xs))} et ^ -^{3,{xp.Xs)\y ~ gixp.Xs))} (4.9) 
oxp dxs 

où Jp et Jg sont les matrices jacobiennes de la fonction g par rapport aux variables Xp et Xs '■ 

n ( ^^ d{9{Xp,Xs))i . (T ( ^^ d{9{Xp,Xs))i 

(Jp(xp,xs))»,j - — mSî^j — (Js(xp,xs))*j - — dîSïsïi — 

Calcul de la matrice jacobienne Jp 

Pour calculer la matrice jacobienne Jp, on effectue un développement de Taylor à l'ordre 1 de la fonction g par rapport 
à Xp- La matrice recherchée est telle que g(xp + ^Xp^ Xs) = ^(Xp? Xs) + Jp<^Xp + ^(ll'^XplP)- Of; on a : 

5(Xp + ^Xp, Xs) = Kapt - EiAo-^E* (B'^)"i((B'=)-i + Xp + <5Xp + X,)-i(Xp + ^Xp + X,)F° (4.11) 

• On commence par effectuer un développement de ((B'^)^^ + Xp + JXp + Xg)^^ à l'ordre 1 en utilisant le lemme 
d'inversion matricielle^. On obtient : 

{{W)-^ + Xp + (5Xp + Xg)-i = [I - <5XpM-i] + 0(||<5Xp||2) (4.12) 

• On revient alors à l'expression de g{xp + ^Xpj Xs) pour obtenir le développement de Taylor à l'ordre 1 de la fonction 
g par rapport à Xp- En utilisant l'expression de Xp en fonction de Xp (voir Partie 2.1.6, page 17), on obtient : 

aiXp + ^Xp. Xs) = ff(Xp, Xs) - EiAq-^E* (B'=)-iM-iHPDiag{(5xp}GP(I - M-i(Xp + Xg))^" + 0(||<5xpll') (4.13) 

• On en déduit l'expression de la matrice jacobienne Jp : 

Jp = -EiAq ^E* (B'=)-iM-iHPDiag {GP(I - M-i(Xp + Xg))V^°} (4.14) 
(en effet, pour deux vecteurs W\ et Wi de même taille, on a : Diag{î«i}ti;2 = Diag{it;2}wi) 
Calcul de la matrice jacobienne Jg 

En reprenant la même démarche que pour le calcul de Jp, on obtient l'expression de la matrice jacobienne Jg : 

3 

Jg = -EiAô^E* (B'=)-iM-i Y. HiDiag {G?(I - M-i(Xp + Xg))^^} (4.15) 
Retour à l'expression du gradient 



2. Lemme d'inversion matricielle : (P + QRS)"! = P^^ - P"iQ(R"i + SP-iQ)-iSP-l 
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Les gradients selon Xp et Xs sont égaux à : 

1^ = {Diag {GP(I-M-i(Xp + X,))yo} {U^r{M-^)H{Bn-'VE,{A^^)^E\{y - g(xp, Xs))} (4.16) 



— = 3? 1^ Diag {g? (I - M-i(Xp + X,))yo} (Hr)*(M-i)t((B'^)-i)tE2(A„-i)tE* (y - g{Xp, X.)) | (4.17) 

Lorsque l'on considère plusieurs fréquences uj et plusieurs positions de la source k, le gradient correspond à la somme 
des expressions obtenues pour chaque couple {w, fc}. On a donc : 

Diag {gp (I - M^-i(Xp + X,))V-,%} (H p)« 
Eli Diag {Gf(I-M^-i(Xp + X,))y^" ,} (Hf )*_ 

(Mji)t((B^)-i)tE2((A^,p,,)ôi)tEl[y„,,-5^,,(Xp,X.)]| (4.18) 

4.4 Analyse de la méthode d'optimisation pixel par pixel 

La méthode d'optimisation pixel par pixel consiste à faire décroître la valeur du critère en ne faisant évoluer qu'une 
seule composante de l'espace d'état à chaque étape de minimisation : pour un point initial donné, on considère un pixel de 
la zone d'étude et l'on fait évoluer une de ses caractéristiques (xp ou Xs)- A priori, nous cherchons à minimiser le critère 
selon l'axe considéré. Après optimisation de la caractéristique, le contraste obtenu fait office de nouveau point initial et 
l'on passe à une autre composante. 

Pour certains problèmes d'imagerie à ondes diffractées, il est possible de définir de telles variations de façon analytique 
[19] ce qui permet de déterminer un minimiseur exact de façon rapide. Cependant, dans le cadre de notre étude, le fait que 
nous travaillons en multifréquentiel et que nous recherchons deux caractéristiques en chaque pixel rend une telle démarche 
inenvisageable. Pour chaque composante considérée, nous rechercherons donc un pas de progression vérifiant les conditions 
de Wolfe. 

Nous détaillons ci-dessous les étapes de calcul associées à cette méthode d'optimisation (calcul de la valeur du critère et 
du gradient au point initial et après variation en un pixel) ainsi que les coiits de calculs correspondants. Nous comparons 
ensuite cette méthode avec les méthodes de type gradient. 

4.4.1 Calculs à un point initial 
Calcul du critère 



C{Xp,Xs) = 2 EE 11^"''^ 

fe 

Lil fe 

avec M„ = (B^)-i+Xp + Xs 

On remarque que pour calculer le critère à un point initial, deux systèmes linéaires doivent être résolus pour chaque 
fréquence et chaque position de la source. Or, les matrices normales de ces systèmes linéaires sont {A^^p^s)o et M^^. Elles 
ne dépendent donc pas de la position de la source et la première est indépendante du contraste. Afin de gagner en coût 
de calcul, nous choisissons d'exploiter le fait qu'une même matrice normale est commune à plusieurs systèmes linéaires en 
procédant de la façon suivante : 

Lors d'une phase d'initialisation de la méthode : 



dC 

dx 



ac 



-5i^,fe(Xp,Xs)|P 

-«fe)capt + Ei(A„,p,,)o-iE*(B^)-iMji(Xp + X,)l/^%||2 (4.19) 
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- pour chaque fréquence, on effectue la factorisation LU de la matrice (A„_p_s)o : (At^_p_s)o = p s)oU(a^ p s)o 

(coût de calcul : 0{n\)) ; 

Pour chaque point initial considéré : 

- pour chaque fréquence, on effectue la factorisation LU de la matrice M;^ : M^^ = Lm„Um„ (coût de calcul : 0{n%[)) ; 

- pour chaque fréquence et chaque position de la source, on résout quatre systèmes linéaires triangulaires (coût de 
calcul : 2 x O(n^) + 2 x 0{nlj)). On pose : Wo = V^^L:^JXp + ^sW^^- 

{uA (rcsp. um) désigne la taille de la matrice (At^^p^s)o (resp. M^j).) 

Calcul du gradient 



Diag{GP(V:,%-«;o}(HP)* 
E?=i Diag {g?{V°^,-wo} (H?)*_ 

(Mji)t((B^r')^E2((A„,p,,)ôi)tEl[î/<,,fc -5„,,(Xp,x.)]| (4.20) 

Le calcul du gradient du critère à un point initial nécessite la résolution de deux systèmes hnéaires supplémentaires 
pour chaque fréquence et chaque position de la source. Or, les matrices normales sont {A^^p^g)^ et Mjj^. Nous exploitons 
donc à nouveau la décomposition LU des matrices (A^. p _,)o et pour se ramener à la résolution de quatre systèmes 
linéaires triangulaires pour chaque fréquence et chaque position de la source. 



dC 
dx 



dC 
9xs 



ui k 



4.4.2 Variation d'une caractéristique en un pixel 



Calcul du critère 

Si l'on considère un pixel et que l'on fait varier la caractéristique Xp de ce pixel d'un pas a, la matrice de contraste Xp 
devient Xp + ah^'gi où g^ est un vecteur ligne et h^* est un vecteur colonne. En utilisant le lemme d'inversion matricielle, 
on obtient : 



9u.,k{Xp + SXp,Xs) = g.,k{Xp,Xs) - Ei(A„,p,,)o-iE*(B^)-iM-ihP (^^ + g^M-^h^) gP(y^% - wo) (4.21) 

~- V " 

=Ai(a), un scalaire 

Par conséquent, la résolution de deux systèmes linéaires liée au choix du pixel (et donc au choix des vecteurs gj et 
h^') mais indépendante de la valeur du pas a est nécessaire pour chaque fréquence. Les matrices normales étant (A^^ p ,j)o 
et Mjj, on utilise à nouveau la factorisation LU de ces deux matrices et on résout quatre systèmes linéaires triangulaires. 
Un changement de la valeur de a n'implique pas la résolution de systèmes linéaires supplémentaires puisque ce coefficient 
n'influe que sur la valeur de Ai (a). On pose : wi = L^^ h^*. 

Le raisonnement est similaire pour Xs mais le nombre de systèmes linéaires triangulaires à résoudre est multiplié par 
trois. 

Calcul du gradient 

En utilisant le lemme d'inversion matricielle, on obtient : 



Diag {gp(1/^', -wo- «;iAi(a)gP(V;% - u;o))} (H?)* 
Eti Diag {G?(y,% -wo- îOiAi(a)gP(K,% - wo))} (Hf)*_ 

(Mt,)-i(/-îi;iAi(a)gP)t((BS)-i)tE2((A„,p,«)ô^)tE*[y„,,-5„,fe(Xp + ^Xp,X«)]| (4-22) 

Par conséquent, la résolution de deux systèmes linéaires qui dépendent cette fois-ci de la valeur du pas a est nécessaire 
pour chaque fréquence et chaque position de la source. Les matrices normales étant (A^,^p^s)q et Mj^, on utilise à nouveau 
les facteurs L(A„_p_,)o, U(A„,p,^)o, Lm„ et Um„ correspondants. 

Le raisonnement est similaire pour Xs- 



VC(xp + ^Xp,X.) = EE^ 

oj k 
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4.4.3 Comparaison avec les méthodes de type gradient 

On résume tout d'abord dans le tableau ci-dessous les étapes de calcul les plus coûteuses associées à la méthode 
d'optimisation pixel par pixel. On ne tient pas compte des calculs effectués lors de la phase d'initialisation de la méthode 
(décomposition LU de la matrice {Ai^^p^s)o) car ils n'entrent pas dans le processus itératif. 



Evolution de Xp un pixel 


Coût 


Pour le point initial 


Décomposition LU des matrices M^^ 

Nf décompositions LU 
Calcul du critère au point initial 

ANfNk résolutions de syst. lin. triangulaires 
Calcul du gradient au point initial 

ANfNk résolutions de syst. lin. triangulaires 
Choix du pixel 

ANf résolutions de syst. lin. triangulaires 


Nf X Oinlj) 
Nf xNk{2x 0{n\) + 2 x 0{nlj)) 
Nf xNk{2x 0{n\) + 2 x 0{nl,)) 
Nf (2 X 0{n\) + 2 X 0{nl,)) 


Pour chaque pas testé 


Calcul du gradient 

ANfNk résolutions de syst. lin. triangulaires 


Nf xNk{2x 0{n\) + 2 x 0{nlj)) 




Evolution de Xs en un pixel 


Coût 


Pour le point initial 


Décomposition LU des matrices M^^ 

Nf décompositions LU 
Calcul du critère au point initial 

ANfNk résolutions de syst. lin. triangulaires 
Calcul du gradient au point initial 

ANfNk résolutions de syst. lin. triangulaires 
Choix du pixel 

12Nf résolutions de syst. lin. triangulaires 


Nf X 0{nl,) 
Nf xNk{2x 0{n\) + 2 x 0{nl,)) 
Nf xNk{2x 0{n\) + 2 x 0(n|^)) 
SxNf{2x 0(n\) + 2 X 0(n|^)) 


Pour chaque pas testé 


Calcul du gradient 

ANfNk résolutions de syst. lin. triangulaires 


Nf xNk{2x 0{n\) + 2 x 0{nlt)) 



Comme pour la méthode d'optimisation pixel par pixel, une méthode de type gradient nécessite le calcul de la valeur 
du critère et du gradient en différents points de l'espace de représentation. Pour chaque point considéré, on est alors 
amené à résoudre plusieurs systèmes linéaires dont les matrices normales sont (A„^p_s)o (ou {Au,p^s)o) et M^^ (ou Mj,). 
On procède donc de la façon suivante : 

Lors d'une phase d'initialisation de la méthode : 

- pour chaque fréquence, on effectue la factorisation LU de la matrice {A.i^^p^g)o (coût de calcul ; Nf x 0{n\)) ; 
Pour chaque point de l'espace de représentation considéré : 

- pour chaque fréquence, on effectue la factorisation LU de la matrice M^^ (coût de calcul : Nf x ©(n^)) ; 

- pour calculer la valeur du critère, on utilise les facteurs LU pour se ramener à la résolution de ANfNk systèmes 
linéaires triangulaires (coût de calcul : N f x Nk{2 x 0{n\) + 2 x 0(n|^))) ; 

- pour calculer le gradient du critère, on utilise les facteurs LU pour se ramener à la résolution de ANfNk systèmes 
linéaires triangulaires supplémentaires (coût de calcul : Nf x Nk{2 x O(n^) + 2 x 0{n%[))). 

Etant donne que l'algorithme de Moré et Thuente ne nécessite que peu d'essais pour obtenir un pas de progression 
vérifiant les conditions de Wolfe, le temps gagné à chaque itération en utilisant la méthode d'optimisation pixel par pixel 
n'est pas significatif par rapport à une méthode de type gradient. 

De plus, pour les méthodes de type gradient, la direction de recherche est choisie en fonction des variations locales du 
critère dans l'espace de représentation, ce qui n'est pas le cas pour la méthode d'optimisation pixel par pixel (à chaque 
itération, on se restreint aux variations du critère le long d'un des axes principaux). On s'attend donc à ce qu'une méthode 
de type gradient converge en un nombre d'itérations beaucoup plus réduit. 

L'utilisation d'une méthode de type gradient semble donc plus adaptée à notre problème. 

4.5 Choix d'une méthode pour définir la direction de recherche 

Pour définir une direction de recherche, nous avons proposé : 

- le gradient conjugué non linéaire ; 
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- l'algorithme L-BFGS ; 

- une méthode d'optimisation pixel par pixel. 

Une analyse de la méthode d'optimisation pixel par pixel nous a montré qu'elle serait moins efficace que les autres pour 
le problème traité ici. Pour faire un choix parmi les méthodes restantes, nous avons appliqué les algorithmes d'inversion 
pour chacune de ces méthodes (nous avons utilisé le milieu test de petite taille présenté Partie 2.5 page 20 ; les données 
mesurées ne sont pas bruitées et il n'y a pas de terme de régularisation). Nous présentons sur la Figure 4.1 l'évolution du 
critère pour chaque cas au cours des 2000 premières itérations ; nous y incluons également l'évolution du critère obtenue 
avec la méthode de plus forte pente. 

Remarque : Il existe plusieurs variantes de l'algorithme du gradient conjugué non linéaire. Nous avons retenu l'algo- 
rithme de Polak-Ribière qui s'avère généralement plus efficace [4]. 



Steepest Descent 

Gradient Conjugué - Polak-Ribière 
L-BFGS - m = 3 
L-BFGS - m = 7 




Figure 4.1 - Choix d'une méthode pour définir la direction de recherche. Les courbes montrent l'évolution temporelle du 
critère pour quatre méthodes de type gradient : plus forte pente, gradient conjugué non linéaire (Polak Ribière), L-BFGS 
à l'ordre 3 et à l'ordre 7. 



Hormis l'algorithme de plus forte pente, les évolutions temporelles de la valeur du critère obtenues avec les différentes 
méthodes sont similaires. Etant donné que l'algorithme L-BFGS nécessite de stocker en mémoire un nombre de variables 
plus important que l'algorithme du gradient conjugué (l'ordre m désigne le nombre de directions de recherche utilisées 
pour définir la direction de recherche à l'itération courante) , nous choisissons finalement de retenir l'algorithme du gradient 
conjugué non linéaire pour définir la direction de descente à chaque itération. 



4.6 Premiers résultats 



Nous avons appliqué l'algorithme d'inversion au milieu de petite taille présenté dans la Partie 2.5, page 20. Les 
données mesurées correspondent aux données obtenues par résolution du problème direct auxquelles nous avons ajouté 
un bruit blanc gaussien tel que le rapport signal à bruit soit égal à 30dB. Le critère comprend désormais deux termes de 
régularisation : 

- un premier terme de rappel aux caractéristiques de la terre (rappel à zéro), on utilise la norme Ll. On pénalise ainsi 
les valeurs de contraste d'amplitude trop grande. Le coefficient associé à ce terme est fixé à 10" 

- un second terme de différence entre pixels voisins, on utilise une norme L1L2 {Clil2 = •v/llxP + où le vecteur x 
est la concaténation des vecteurs Xp et Xs)- Le paramètre S est fixé à 10^ et le coefficient associé à ce terme est fixé 
à 10-18. 
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Nous présentons sur la Figure 4.2 les résultats obtenus pour deux initialisations différentes : pour la première initia- 
lisation, les caractéristiques de la zone d'étude sont égales à celles de la terre (les deux contrastes Xp 6t Xs initiaux sont 
nuls dans toute la zone d'étude) et pour la seconde, les contrastes sont égaux à la solution recherchée. Nous montrons les 
cartes obtenues ainsi que l'évolution temporelle de la valeur du critère (critère total et terme d'adéquations aux données) 
et de la norme du gradient pour ces deux initialisations. 

On remarque tout d'abord un problème de convergence : les résultats présentés ici ont été obtenus après 30000 itérations 
(soit 16 heures de calcul environ) et la convergence n'a toujours pas été atteinte. Pour l'initialisation de la zone d'étude à 
des contrastes nuls (initialisation aux caractéristiques de la terre), le critère prend des valeurs de plus en plus proches de 
celles obtenues en initialisant à la solution mais il diminue de façon très lente. 

On note également que pour l'initialisation à des contrastes nuls, les valeurs de contraste obtenues après 30000 itérations 
sont élevées mais encore très éloignées de celles du béton : le contraste Xp atteint la valeur maximale de 9,9.10^m^/s^ 
au lieu de 15, 9.10^m^/s^ et le contraste Xs atteint la valeur maximale de 6,6.10^m^/s^ au lieu de 4, S.lO^m^/s^. Les 
contrastes obtenus en partant de contrastes nuls sont donc encore loin de la solution recherchée. 

Si l'on s'intéresse aux contrastes obtenus en initialisant les contrastes à zéro au bout de 10000 itérations (nombre d'ité- 
rations à partir duquel le critère évolue lentement) , on remarque que les cartes obtenues sont déjà proches de celles obtenues 
au bout de 30000 itérations et que les valeurs de contrastes atteintes sont du même ordre de grandeur : 6,6.10^m^/s^ 
pour Xp et 4, 2.10^m^/s^ pour Xs- La lenteur de l'évolution du critère semble donc coïncider avec la forte amplitude des 
contrastes obtenus. 

4.7 Accélération de la convergence avec un changement de variable 

4.7.1 Expression du critère en fonction d'un autre jeu de variables caractéristiques 

Les observations précédentes semblent mettre en évidence un problème de sensibilité du critère vis-à-vis des variations 
de contraste Xp et Xs lorsque ces derniers prennent des valeurs élevées. En effet, pour l'initialisation aux caractéristiques 
de la terre (contrastes nuls), on observe une rapide décroissance du critère lors des premières itérations puis le critère 
décroît plus lentement et les valeurs des contrastes correspondants sont plus élevées, bien qu'encore éloignées des valeurs 
recherchées. De même, pour l'initialisation à la solution, Xp et Xs prennent des valeurs élevées et le critère évolue lentement 
dès les premières itérations. 

Afin d'améliorer la sensibilité du critère, nous proposons d'effectuer un changement de variable : au lieu d'exprimer 
les termes d'adéquation aux données et de régularisation du critère en fonction des variables Xp et Xs, on les exprime 
en fonction des variables Cp et as qui sont choisies de sorte que de faibles variations de ap et CTs induisent de fortes 
variations de Xp et Xs pour des valeurs de contraste élevées. Ainsi, le critère à minimiser reste identique mais les directions 
de recherches sélectionnées diffèrent d'une variable à l'autre (le changement de variable a une incidence sur le calcul du 
gradient). 

Les différentes variables utilisées ont été présentées dans la Partie 2.4, page 19. Nous avons effectué ces changements 
de variable en reprenant le même cas d'étude que précédemment : 

- on considère le milieu de petite taille ; 

- les données mesurées sont obtenues par résolution du problème direct avec ajout de bruit blanc gaussien (le rapport 
signal à bruit est égal à 30dB) ; 

- le critère comprend deux termes de régularisation : un premier terme de rappel aux caractéristiques de la terre 
utilisant la norme Ll avec un coefficient de pondération de 10""'^^ et un second terme de différence entre pixels 
voisins utilisant une norme L1L2 dont le paramètre S est fixé à 10^ avec un coefficient de pondération de 10~^^. 

On représente sur la Figure 4.3 les évolutions du critère obtenues pour les différents changements de variable lors 
des 5000 premières itérations. On remarque tout d'abord que pour les différents changements de variable proposés, le 
critère décroît plus rapidement qu'avec les variables Xp et Xs- Cela confirme le fait qu'un problème de sensibilité du critère 
vis-à-vis de ces variables ralentissait la convergence de l'algorithme dans le cas précédent. La décroissance la plus rapide 
est observée avec l'utilisation des variables ap — Invp et CTs — Invg. Nous conserverons donc ce changement de variable 
par la suite. 

Remarque sur la prise en compte des contraintes de positivité sur les vitesses : 

Les vitesses de propagation des ondes en pression (vp) et en cisaillement (vg) sont des grandeurs positives. Or, pour 
certains changements de variable, on risque de passer par des valeurs négatives de Vp et Vg, ce qui peut engendrer un 
comportement pathologique de l'algorithme. C'est le cas par exemple lors de l'utilisation des variables Up = ^ et CTs = ^ 
qui ne devraient pas prendre de valeurs négatives. 



D. Vautrin, M. Voorons, J. Hier, Y. Goussard 



61 



Chapitre 4. Méthodes inverses fondées sur la formulation primale du problème direct 




j ^ ^ ^ ^ 1 ,0-1 ^ ^ ^ ^ ^ 

8 10 12 14 16 18 2 4 6 8 10 

Temps (h) Temps (h) 



(c) Evolution temporelle du critère (critère total et terme d'adé- (d) Evolution temporelle de la norme du gradient 

quation aux données) 



Figure 4.2 - Premiers résultats obtenus sur le milieu de petite taille 
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10"" I ' ' 1 1 ' ' 1 ' 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 

Temps (h) 

Figure 4.3 - Evolution temporelle du critère pour les différents changements de variable proposés en initialisant aux 
caractéristiques de la terre (en noir : cTp ^ — Xp,s, en bleu : s ~ Vp,s, en vert : Cp = en rouge : ap^g — In î'p.s) 



Dans ce cas, il est possible de respecter les contraintes de positivité tout en évitant une incidence notable sur le 
comportement de l'algorithme d'inversion en effectuant un changement de variable supplémentaire. On propose d'utiliser 
les variables ^pp ou ips telles que : 

Ainsi, on s'assure que la contrainte de positivité est respectée et si le coefficient e est choisi suffisamment petit, le chan- 
gement de variable n'a qu'une faible incidence au-delà d'un certain seuil strictement positif (la fonction est proche de la 
fonction identité). On représente sur la figure 4.4 le tracé de la fonction utilisée pour passer de Tpp^s à (Jp^s- 



5 



f(x) = 1/2[{ + + X] 

f(x) - X 











-1 ' • • • ' • • • • 1 

-5 -4 -3 -2 -1 1 2 3 4 5 



X 

Figure 4.4 - Tracé de la fonction utilisée pour vérifier la contrainte de positivité sur les vitesses avec e — 1 (trait plein) 
et comparaison avec la fonction identité (pointillés) 

Pour le changement de variable retenu, les variables Cp et a g sont forcément associées à des valeurs de Vp et Vg positives 
(on a Vp — expcTp et î;^ — expcr^). Un tel changement de variable est donc inutile dans ce cas. 
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Remarque sur l'affichage des résultats : 

Sur la Figure 4.2, nous avons présenté les premiers résultats obtenus sur le milieu de petite taille en affichant les 
contrastes Xp 6t Xs- Nous afficherons maintenant les cartes correspondant aux nouvelles variables utilisées, c'est-à-dire 
tTp = In Vp et (7 s = In «s . 

4.7.2 Pénalisation d'un autre jeu de variables 

Pour les premiers résultats obtenus, les deux termes de régularisation du critère portaient sur les contrastes Xp st Xs- 
Cependant, cela peut induire des problèmes de conditionnement (on a Xp.s = Gxp2ap^s — exp2crp_s,o pour ap^s = hiVp,^). 
Par exemple, des variations de <Tp ou CTs en un pixel autour d'une grande valeur auront une forte incidence sur le terme 
de rappel aux caractéristiques de la terre contrairement à des variations autour d'une valeur plus faible. C'est pourquoi il 
peut être préférable d'appliquer la régularisation à un autre jeu de variables. 

Nous avons lancé l'algorithme d'inversion sur le milieu de petite taille en portant la régularisation sur les variables 
CTp = In Vp et as = In Vs et en exprimant le critère en fonction de ces mêmes variables. Comme précédemment, nous utilisons 
les données obtenues par résolution du problème direct auxquelles nous avons ajouté un bruit blanc gaussien (rapport 
signal à bruit égal à 30dB) et le critère comprend deux termes de régularisation : 

- un terme de rappel aux caractéristiques de la terre pour lequel on utilise la norme Lf avec un coefficient de pondé- 
ration égal à 10"^'^ ; 

- un terme de différence entre pixels voisins pour lequel on utilise une norme Lf L2 avec un paramètre S égal à 0, f et 
un coefficient de pondération égal à fO~^^. 

Nous considérons que l'algorithme est arrivé à convergence lorsque la norme du gradient divisée par le nombre d'incon- 
nues (égal à 300 ici) est inférieure à f O^"^** pendant 50 itérations successives. Nous présentons sur la Figure 4.5 les résultats 
obtenus en initialisant les caractéristiques de la zone d'étude à celles de la terre d'une part et à la solution recherchée 
d'autre part. Nous montrons les cartes obtenues ainsi que l'évolution temporelle de la valeur du critère et de la norme du 
gradient. 

Le fait d'exprimer le critère en fonction de ap = In Vp et CTs = In Vs et d'utiliser ces variables pour la régularisation 
a permis d'améliorer le comportement de l'algorithme d'inversion. La sensibilité du critère vis-à-vis des variables opti- 
misées est améliorée puisqu'en moins de six heures, le critère a maintenant convergé vers la même valeur pour les deux 
initialisations. 

On note également que pour les deux initialisations, les deux cartes obtenues sont similaires. Cependant, les valeurs 
maximales n'atteignent pas les valeurs recherchées (7, 6 au lieu de 8, 3 pour In Vp et 7, 4 au lieu de 7, 7 pour In ii^). Cela est 
dii aux termes de régularisation du critère qui tendent à diminuer les amplitudes (on remarque que cette diminution des 
valeurs maximales n'affecte quasiment pas la valeur du terme d'adéquation aux données du critère pour l'initialisation à 
la solution). 

4.8 Introduction des fréquences de façon progressive 

Jusqu'à maintenant, nous avons utilisé l'algorithme d'inversion en introduisant dès le départ toute l'information fré- 
quentielle. Or, il peut être préférable d'introduire les informations contenues dans les basses fréquences dans un premier 
temps puis d'ajouter les informations contenues dans les plus hautes fréquences de manière progressive : les basses fré- 
quences apportent des informations sur les variations spatiales lentes et permettent d'obtenir une image lisse du milieu ; 
l'introduction des plus hautes fréquences permet ensuite d'affiner cette image [9]. 

Nous avons testé cette démarche en procédant de la façon suivante : tout d'abord, nous n'avons introduit que les 
informations associées à la plus basse fréquence (46,7 Hz). L'algorithme étant arrivé à convergence, nous avons ensuite 
ajouté les informations associées à la fréquence suivante (93,3 Hz) puis nous avons procédé de la même manière pour les 
fréquences plus élevées jusqu'à prendre en compte toutes les fréquences (f5 fréquences allant de 46,7 Hz à 700 Hz). 

Pour chaque groupe de fréquences, nous considérons que l'algorithme est arrivé à convergence lorsque la norme du 
gradient divisée par le nombre d'inconnues est inférieure à 5.f0~^'^ pendant 50 itérations successives. Lorsque ce critère 
d'arrêt est vérifié, on utilise les cartes obtenues pour initialiser l'algorithme d'inversion que l'on relance en introduisant 
une fréquence supplémentaire. 

Les valeurs des coefficients intervenant dans les deux termes de régularisation sont les mêmes que précédemment (le 
coefficient de pondération associé au terme de rappel aux caractéristiques de la terre est égal à 10"^'^ et le coefficient 
de pondération associé au terme de différence entre pixels voisins est égal à fO""'^^ ; le paramètre 6 est égal à 0, f). Les 
résultats obtenus sont présentés sur la Figure 4.6. Nous affichons également sur la Figure 4.7 les données mesurées par 
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(a) Cartes obtenues en initialisant aux caractéristiques de la terre (à gauche : 
InVp, à droite : In «s ) 
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(b) Cartes obtenues en initialisant à la solution (à gauche : In-ï^p, à droite 
In Vs) 



Initialisation : uniforme, terre 
Initialisation : uniforme, terre (AD) " 
- Initialisation : solution 
■ Initialisation : solution (AD) 



- Initialisation : uniforme, terre 
■ Initialisation : solution 




3 

Temps (h) 



(c) Evolution temporelle du critère (critère total et terme d'adé- 
quation aux données) 



Temps (h) 

(d) Evolution temporelle de la norme du gradient 



Figure 4.5 - Résultats obtenus sur le milieu de petite taille en utilisant les variables ap = In Vp et a s — In Vs 
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les capteurs lorsque les cartes correspondent à la solution recherchée (sans et avec ajout de bruit blanc) et aux résultats 
obtenus pour les deux initialisations. 

Les résultats obtenus avec l'introduction progressive des fréquences restent satisfaisants : l'algorithme arrive à conver- 
gence et les cartes obtenues sont similaires à celles obtenues précédemment (voir Figure 4.5). De plus, la comparaison des 
données reconstruites montrent que pour les deux initialisations, les mesures sont semblables. On remarque cependant que 
l'on ne parvient pas à retrouver le contenu hautes fréquences des données utilisées pour l'inversion. Cela s'explique par le 
fait que les composantes hautes fréquences sont de faible amplitude et n'ont donc qu'une faible influence sur la valeur du 
critère. 

L'intérêt principal de cette démarche est la diminution du temps de calcul : sans l'introduction progressive des fré- 
quences, il fallait environ trois heures de calcul pour l'initialisation aux caractéristiques de la terre et presque six heures 
de calcul pour l'initialisation à la solution. Il faut maintenant une heure et demie environ à l'algorithme pour arriver à 
convergence pour les deux initialisations. Cela s'explique surtout par le fait que le coiit de calcul par itération pour les 
premiers groupes de fréquences est plus faible (le temps de calcul du critère et du gradient est quasiment proportionnel 
au nombre de fréquences considérées). 

4.9 Résultats obtenus sur le milieu de taille intermédiaire 

Nous présentons maintenant les résultats obtenus sur le milieu de taille intermédiaire présenté dans la Partie 2.5, page 
20 en tenant compte des conclusions auxquelles nous sommes parvenus d'après les résultats obtenus sur le petit milieu : 

- le critère à minimiser s'exprime en fonction des variables CTp — In Vp et as — Invs', 

- les termes de régularisation du critère portent sur cette même variable ; 

- nous introduisons les fréquences de manière progressives, des basses fréquences vers les hautes fréquences. 
Comme pour le petit milieu, les données mesurées sont obtenues par résolution du problème direct puis ajout de bruit 

blanc gaussien (le rapport signal à bruit est égal à 30dB) et le critère comprend deux termes de régularisation : un terme 
de rappel aux caractéristiques de la terre pour lequel on utilise la norme Ll avec un coefficient de pondération égal à 10~^^ 
et un terme de différence entre pixels voisins pour lequel on utilise une norme L1L2 avec un paramètre S est fixé à 0, 1 et 
un coefficient de pondération égal à 10^^". 

Nous avons considéré que l'algorithme était arrivé à convergence lorsque la norme du gradient divisée par le nombre 
d'inconnues (environ 1400 ici) est inférieure à 5.10^^"^ pendant 50 itérations successives. 

Nous présentons sur la Figure 4.8 les résultats obtenus pour deux initialisations différentes (pour la première initialisa- 
tion, les caractéristiques de la zone d'étude sont égales à celles de la terre et pour la seconde, les contrastes sont égaux à la 
solution recherchée). Nous présentons également sur la Figure 4.9 une comparaison des données mesurées par les capteurs 
lorsque les cartes correspondent à la solution recherchée (sans et avec ajout de bruit blanc) et aux résultats obtenus pour 
les deux initialisations (les sismogrammes sont maintenant représentés en niveaux de gris étant donné le nombre élevé de 
capteurs) . 

On remarque tout d'abord que l'algorithme arrive bien à convergence : pour les deux initialisations, les cartes obtenues 
sont similaires et le critère converge vers la même valeur. Les cartes obtenues montrent que la forme de l'objet diffractant 
est plutôt bien reconstruite. On remarque cependant que les contours sont plus marqués sur les cartes des a s = In î^s- Sur 
les différentes cartes, les deux régions (terre et béton) ne sont pas parfaitement homogènes mais les valeurs caractéristiques 
correspondantes sont proches de celles du milieu recherché. Comme pour le milieu de petite taille, la comparaison des 
données reconstruites montrent que les mesures obtenues sont semblables pour les deux initialisations et que l'on retrouve 
bien les données utilisées pour l'inversion. 

Il a fallu davantage de temps (presque deux jours) par rapport au milieu de petite taille pour arriver à convergence 
mais plusieurs modifications pourraient encore être apportées afin de réduire le temps de calcul comme par exemple : 

- pouvoir se ramener à des données non redondantes afin de réduire la quantité d'informations traitées (éliminer 
certaines fréquences ou certains capteurs) ; 

- modifier l'algorithme de minimisation en ayant recours au préconditionnement ou à des décompositions LU incom- 
plètes afin d'arriver plus rapidement à convergence pour chaque groupe de fréquences considéré. 
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(a) Cartes obtenues en initialisant aux caractéristiques de la terre (à gauche : 
InVp, à droite : InVg) 
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(b) Cartes obtenues en initialisant à la solution (à gauche : InîJp, à droite 
In Vs ) 



Inilialisation : unilorme, terre 
initialisation : Uniterme, terre (AD) 
initialisation : soiution 
initialisation : soiution (AD) 




- initialisation : uniforme, terre 
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(c) Evolution temporelle du critère (critère total et terme d'adé- 
quation aux données) 



(d) Evolution temporelle de la norme du gradient 



Figure 4.6 - Résultats obtenus sur le milieu de petite taille en utilisant les variables ap = In Vp et <7s = In Vs et en 
introduisant les fréquences de façon progressive 
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(a) Pour chaque fréquence considérée, amplitude des données bruitées utilisées pour l'inversion (rouge pointillé) et 
amplitude des mesures correspondant à la solution recherchée (rouge), aux cartes obtenues en initialisant à la terre (vert) 
et aux cartes obtenues en initialisant à la solution (bleu) en chaque capteur 
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(b) Sismogrammes reconstruits à partir des cartes obtenues (c) Sismogrammes reconstruits à partir des cartes obtenues 
en initialisant aux caractéristiques de la terre en initialisant à la solution 
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(d) Sismogrammes correspondant à la solution recherchée 



Figure 4.7 - Comparaison des données capteurs correspondant aux cartes solutions (sans et avec ajout de bruit blanc) 
et aux cartes obtenues après inversion (initialisation à la terre et à la solution) dans le cas où la source est positionnée 
à gauche de l'objet diffractant (première position sur la Figure 2.2, page 21). On affiche les données dans le domaine 
fréquentiel et dans le domaine temporel (sismogrammes). 
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(b) Cartes obtenues en initialisant aux caractéristiques de la terre (à gauche : InVp, à droite 
Invs) 



(c) Cartes obtenues en initialisant à la solution (à gauche : Infp, à droite ; In «a) 
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(d) Evolution temporelle du critère (critère total et terme d'adé- 
quation aux données) 



(e) Evolution temporelle de la norme du gradient 



Figure 4.8 - Résultats obtenus sur le milieu de taille intermédiaire 
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(a) Pour chaque fréquence considérée, amplitude des données bruitées utilisées pour l'inversion (rouge pointillé) et 
amplitude des mesures correspondant à la solution recherchée (rouge), aux cartes obtenues en initialisant à la terre 
(vert) et aux cartes obtenues en initialisant à la solution (bleu) en chaque capteur 
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(b) Sismograrnmcs reconstruits à partir des cartes obtenues 
en initialisant aux caractéristiques de la terre 



(c) Sismogramrncs reconstruits à partir des cartes obtcimcs 
en initialisant à la solution 
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(d) Sismogrammes correspondant à la solution recherchée 



- Avec les cartes obtenues en initialisant à la terre 

- Avec les cartes obtenues en initialisant à la solution 

- Avec les cartes correspondant au milieu reciierctié 
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(e) Comparaison des mesures à l'instant t = 25 ms (corres- 
pond au trait noir affiché sur les sismogrammes) 



Figure 4.9 - Comparaison des données capteurs correspondant aux cartes solutions (sans et avec ajout de bruit blanc) 
et aux cartes obtenues après inversion (initialisation à la terre et à la solution) dans le cas où la source est positionnée 
à gauche de l'objet diffractant (septième position sur la Figure 2.3, page 21). On affiche les données dans le domaine 
fréquentiel et dans le domaine temporel (sismogrammes). 



D. Vautrin, M. Voorons, J. Idier, Y. Goussard 



70 



Conclusion et perspectives 



Nous avons abordé l'imagerie de la subsurface en distinguant deux familles de méthodes : la première regroupe les mé- 
thodes faisant intervenir des variables auxiliaires et la seconde comprend les méthodes utilisant une formulation "primale" 

du problème direct. 

Les méthodes faisant intervenir des variables auxiliaires ont l'avantage d'utiliser des formulations pour lesquelles le 
critère et le gradient sont rapides à calculer. Cependant, elles font intervenir un grand nombre de grandeurs à optimiser 
ce qui tend à augmenter le nombre d'itérations nécessaires et donc le temps de calc;ul global. Les méthodes s'appuyant 
sur une formulation primale permettent de minimiser le nombre de variables à optimiser et donc le nombre d'itérations 
nécessaire. Cependant, la formulation utilisée est complexe, ce qui implique un temps de calcul élevé à chaque itération. La 
question du compromis entre nombre de variables à manipuler, nombre total d'itérations et volume de calcul par itération 
est donc critique pour une résolution satisfaisante du problème. 

Nous avons étudié différentes méthodes basées sur une formulation bilinéaire du problème. Toutes font apparaître des 
formes algébriques différentes et ont des qualités qui leur sont propres. Certaines peuvent être appliquées sur une zone 
d'étude, limitant ainsi le nombre d'inconnues, tandis que d'autres bénéficient d'estimateurs simples à calculer. Nous ne 
sommes cependant parvenus à en faire converger aucune, ce qui rend difficile toute comparaison quantitative. Nous pouvons 
tout de même souligner que la méthode du gradient modifie semble être la plus intéressante des formulations bilinéaires, 
puisqu'elle permet d'obtenir les résultats plus proches de la solution après un nombre d'itération donné. La méthode CFSI 
a un important problème de conditionnement croisé de deux de ses estimateurs, tandis que la méthode CSI a un coût de 
calcul par itération plus important que la méthode CM. Finalement, la méthode sans contraste est la plus gourmande en 
place mémoire et semble converger plus lentement que toutes les autres. Dans tous les cas, le critère décroît très lentement 
après quelques centaines d'itérations et aucune méthode semble converger après plusieurs jours de calcul. La lenteur avec 
laquelle le critère décroit peut être expliquée par le mauvais c;onditionnement de la plupart des matrices normales. Ce 
conditionnement se détériore au fur et à mesure des itérations ; il en résulte une diminution du déplacement relatif de la 
solution, ce qui entraîne une stagnation du critère. Les changements de variables proposés ont eu des effets différents, mais 
aucun ne s'est avéré être une solution au problème de convergence de la méthode à laquelle il a été appliqué. Des tests 
supplémentaires ont été menés pour tenter de déterminer la raison des problèmes de convergence rencontrés. Les résultats 
de ces tests sont donnés en annexe de ce document et ne permettent pas de conclure avec certitude sur l'origine de ces 
problèmes de convergence. 

En ce qui concerne la formulation primale, nous avons tout d'abord proposé une formulation optimisée de manière à 
diminuer à la fois le temps et l'espace mémoire nécessaires au calcul du critère et du gradient. Une première analyse des 
méthodes envisagées nous a amené à retenir celle du gradient conjugué non linéaire. Les premiers essais effectués sur le 
milieu de petite taille nous ont conduit à exprimer le critère en fonction d'autres variables que le contraste (In Vp et In Vs) 
et à utiliser ces mêmes variables pour la régularisation. 

Pour les deux familles de méthodes, l'introduction progressive des fréquences nous a permis de diminuer le temps de 
calcul. Dans le cas des méthodes bilinéaires, cette accélération permet seulement d'atteindre plus rapidement le palier de 
convergence lente du critère. Par contre, pour la formulation primale, la procédure d'inversion ainsi obtenue s'est révélée 
plutôt efficace et nous a permis de faire des tests sur le milieu de taille intermédiaire. Plusieurs pistes pourraient cependant 
être explorées pour accélérer davantage l'algorithme d'inversion. 

Au vu des résultats obtenus pour les différentes méthodes abordées, l'utilisation de l'algorithme utilisant la formulation 
primale semble plus adaptée à notre problème. Cependant, des évolutions restent à prévoir. D'une part, l'algorithme doit 
être en mesure de traiter des problèmes de plus grande taille : les milieux traités dans la réalité seront de plus grande 
dimension et avec une résolution spatiale plus fine. D'autre part, il faut l'adapter pour qu'il prenne en compte la présence 
d'une surface libre (pour les problèmes traités jusqu'à maintenant, l'objet diffractant ainsi que la source et les capteurs 
étaient enfouis sous terre). 

Enfin, certains points restent encore en suspens comme, par exemple, la détermination du milieu de référence qui est 
choisi par l'utilisateur et les fréquences nécessaires à une inversion de qualité. 
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Conclusion et perspectives 



Les travaux effectués jusqu'à présent nous ont donc permis de développer une première méthode d'inversion 2D de type 
"cartographie" n'incluant pas d'information a priori concernant les valeurs caractéristiques du béton et la géométrie de 
l'objet recherché. Celle-ci s'est avérée performante sur des jeux de données synthétiques de tailles petite et intermédiaire. 
Dans une certaine mesure, cette méthode constitue donc une « preuve de concept » de la faisabilité de l'inversion de 
données sismiques 2D pour l'examen de fondations de pylônes, et elle servira de base à la suite du projet. Pour aboutir à 
une méthode satisfaisant les objectifs fixés au début de l'étude, les étapes ultérieures seront les suivantes : 

- Poursuite de l'étude de la méthode développée jusqu'ici — L'objectif sera d'en caractériser plus précisément 
les performances, d'en améliorer le comportement numérique et d'élargir la gamme de données qu'elle est capable 
de traiter. 

- Développement de modélisations plus précises du milieu — Il s'agira d'inclure davantage d'informations 
a priori sur la nature de l'objet recherché, notamment la présence de deux zones différentes (la fondation et le 
sous-sol) dans le milieu à imager . Deux approches, réparties entre les institutions impliquées dans l'inversion, sont 
envisagées : l'IRCCyN étudiera les techniques de type contour tandis que l'Ecole Polytechnique de Montréal abordera 
les méthodes de type région. Dans un cas comme dans l'autre, ceci devrait permettre de mieux identifier la forme 
de la fondation et donc contribuer à atteindre les objectifs du projet. 

L'ensemble des méthodes sera testé sur les divers types de données, en fonction de leur disponibilité : données synthétiques 
de plus grande taille avec surface libre, données obtenues avec les maquettes simulées par le LCPC, données réelles recueillies 
lors des campagnes de mesures. 
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Annexe : 
méthode 



Tests sur la convergence de la 
du gradient modifié 



Introduction 

Suite aux résultats obtenus par les méthodes basées sur l'approche bilinéaire, un certain nombre de questions ont été 
posées auxquelles nous donnons des éléments de réponse dans cette annexe. Les trois questions soulevées concernaient : 

- la possibilité de convergence de la méthode GM vers un minimum local ; 

- l'effet de la pénalisation quadratique du champ de vitesse sur le résultat de la reconstruction ; 

- la comparaison entre les variables auxiliaires réelles et estimées. 

Pour obtenir des résultats rapidement nous avons effectué ces tests sur le milieu de petite taille. Les tests ont été effectués 
avec la représentation bilinéaire en utilisant la méthode qui donnait les meilleurs résultats, à savoir la méthode GM avec 
changement de variable en log. 

Test de convergence vers un minimum local 

Pour essayer de déterminer si le critère que nous tentons de minimiser possède des minima locaux, nous avons effectués 
deux tests. Le premier test consiste à observer la variation du critère lorsque l'on fait varier linéairement les valeurs du 
contraste entre la valeur d'initialisation (on initialise à la terre) et l'estimée obtenue en initialisant l'algorithme avec la 
solution. Le deuxième test est identique dans la manière, mais utilise les valeurs du critère obtenues par la méthode 
GM jusqu'à la valeur obtenue lors de l'arrêt de la méthode après 8000 itérations, puis les contrastes sont calculés de 
façon linéaire entre cette estimée et l'estimée obtenue en initialisant l'algorithme avec la solution. Dans les deux cas, les 
variables auxiliaires sont déduites des contrastes en utilisant l'estimateur direct issu de la formulation du GM donnée par 
les équations 3.16. De même, l'expression du critère de la méthode GM que nous calculons est donné par l'équation 3.18. 

Le résultat de ce test est présenté à la figure A. Son interprétation est délicate. Lorsque l'on part de valeurs du contraste 
nulles, ce qui est le cas si on initialise la méthode avec les valeurs caractéristiques de la terre, alors le critère décroit à 
mesure que le contraste se rapproche de la solution. Par contre, si nous partons de l'estimée obtenue après 8000 itérations 
de la méthode GM, alors on peut voir que le critère n'est pas monotone décroissant. Cependant, on ne peut pas conclure 
que la méthode GM décroit obligatoirement vers un minimum local. En effet, la forme du critère est très probablement 
due à un problème de conditionnement de l'estimateur des variables auxiliaires comme le laisse suggérer la courbe (c) 
de la figure A représentant le coût relié aux variables auxiliaires et qui présente la même allure que le critère partant de 
l'estimée de la méthode GM. 

Effet de la régularisation sur la variable auxiliaire 

Les méthodes d'inversion de type bilinéaires que nous proposons dans ce rapport utilisent une fonction de pénalisation 
quadratique sur les variables auxiliaires. Pour mesurer l'effet de cette régularisation sur le résultat, nous avons fait varier 
le paramètre de régularisation associé à cette fonction de pénalisation et nous avons observé les variations résultantes sur 
l'erreur quadratique moyenne du contraste et de la variable auxiliaire, ainsi que sur le temps de calcul. Dans tous les cas 
l'inversion s'arrêtait lorsque l'on atteignait un seuil de la norme du gradient de 2.10^^^. 

Nous avons limité le choix de l'hyperparamètre aux valeurs comprises dans l'intervalle [10"^^; 10"^] car en deçà de 
-|^Q-i8 l'gffgt (jg la pénalisation est négligeable, et au delà de 10~^ on observait une sur-régularisation marquée. 

Ces résultats indiquent que la pénalisation quadratique des variables auxiliaires a un impact significatif sur le nombre 
d'itérations de la méthode, et donc le temps de calcul, ainsi que sur la qualité du résultat de l'estimation de toutes 
les variables : les contrastes en P et S et les variables auxiliaires. Le temps de calcul augmente lorsque la valeur de 
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(c) Coût lié aux variables auxiliaires 

Figure A - Variation de la fonction de coiit de la méthode GM, en faisant varier le contraste linéairement à partir de 
l'initialisation à la terre (a) ou de la valeur estimée avec la méthode GM après 8000 itérations (b). La courbe (c) montre 
la partie du critère liée à l'estimation des variables auxiliaires. 




Valeur de l'hyperparamètre 



Valeur de l'hyperparamètre 



Figure B - Évolution de l'erreur quadratique relative sur les contrastes en P et S en fonction de la valeur de l'hyperpa- 
ramètre 



l'hyperparamètre augmente tandis que l'EQM des différentes variables estimées décroit. Il faut donc faire un compromis 
entre la qualité du résultat souhaitée et le temps de calcul. 



Comparaison des variables auxiliaires réelles et estimées 

Dans la figure E nous présentons des courbes montrant les données du problème direct, les données initiales et les 
données reconstruites par la méthode GM pour différentes fréquences. La figure F compare les cartes des champs de 
vitesse en x et y pour les mêmes fréquences obtenues par la méthode GM. 
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Figure C - Évolution de l'erreur quadratique relative sur la variable auxiliaire en fonction de la valeur de l'hyperparamètre 



Valeur de l'hyperparamètre 



Figure D - Évolution du temps de calcul en fonction de la valeur de l'hyperparamètre 



On constate que les champs de vitesses en x et y pour les basses fréquences sont bien reconstruits tandis que ceux 
correspondant aux hautes fréquences sont mal restitués. 
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Mesures reconstruites - GM (700 Hz), EOM = 4.2859 
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(e) Fréquence 700 Hz 

Figure E - Comparaison entre les données générées par le problème direct et les données reconstruites par la méthode 
GM pour différentes fréquences 
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Estimée en X 





(a) Fréquence 100 Hz 




(b) Fréquence 250 Hz 
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(c) Fréquence 400 Hz 



Figure F - Comparaison entre les cartes du champ de vitesse données par le modèle direct et les cartes reconstruites par 
la méthode GM pour différentes fréquences 
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(o) Fréquence 700 Hz 

Figure F - Comparaison entre les cartes du champ de vitesse données par le modèle direct et les cartes reconstruites par 
la méthode GM pour différentes fréquences (suite)) 
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